Extracting Spectra from NIRCam Grism data¶
Taking IceAge (ERS-1309; PI: McClure) Data as an Example¶
Fengwu Sun (Steward Observatory, University of Arizona, 2022)
(Contact: fengwusun[AT]arizona.edu; Homepage: http://u.arizona.edu/~fengwusun/)
with significant help from Eiichi Egami (UofA), Nor Pirzkal (STScI) and Marcia Rieke (UofA)
(v0 on July 16, 2022)
(v1 on August 10, 2022)
(v2 on October 21, 2022)
Public PDF of APT File (PID 1309): https://www.stsci.edu/jwst/phase2-public/1309.pdf
A Short Overview:¶
In this notebook, we will go through the data reduction and spectral extraction pipeline of the JWST/NIRCam Wide-Field Slitless Spectroscopy (WFSS) mode. If you need basic introduction to this observation mode, I would recommend you to take a look of the JDox Page on NIRCam/WFSS. In this notebook, I assume that the potential users have basic knowledge of NIRCam and NIRCam grism, e.g.:
- There are two independent modules of NIRCam (module A and B);
- Each module has an instantaneous Field of View (FoV) of ~4.4 arcmin$^2$;
- There are two dispersion directions called R and C;
- The dispersion is around 1 nm/pix with a spectral resolution of R~1600 for point sources;
- When it is taking grism spectra at LW (2.4-5.0µm), the SW channel is taking direct imaging data.
Below is the flow chart for grism data analysis:
In this image, all of the items enclosed in red boxes are input fits files you download from MAST, which include (1) direct LW images, (2) direct SW images and (3) LW grism images. Starting from *_rate.fits files (STScI pipeline stage-1 products), you can use direct images in LW and SW to produce final mosaics images and source catalogs (enclosed with green box). These images and source catalogs can be very helpful for your target selections when extracting grism spectra. However, what you actually need to input for grism spectra extraction is (i) astrometric catalog that you want to register SW images onto, and (ii) coordinates (RA and DEC) of targets with interest. As long as you can register the astrometry of the LW grism images, using the spectral tracing models, grism dispersion models and flux calibration functions that I produced, you can extract calibrated 2D grism spectra for any sources you want. We will get to the point on how to do astrometric calibration and spectral calibration later in this notebook.
List of Content¶
1. Import necessary python modules¶
In [24]:
''' Astropy modules (I'm using astropy v5.1) '''
from astropy.stats import sigma_clipped_stats, SigmaClip
from astropy.io import fits, ascii
from astropy.visualization.mpl_normalize import ImageNormalize
from astropy.visualization import MinMaxInterval, PercentileInterval, SqrtStretch, ImageNormalize
from astropy import units as u
from astropy import constants as c
from astropy import wcs, table
from astropy.coordinates import SkyCoord
from astropy.table import Table, Column, vstack
from astropy.time import Time
from astropy.modeling.models import custom_model
from astropy.cosmology import FlatLambdaCDM
cosmo = FlatLambdaCDM(H0 = 70, Om0 = 0.3)
''' Photutils modules (I'm using 1.5.0)'''
from photutils import find_peaks, psf, DAOStarFinder, IRAFStarFinder
from photutils import CircularAperture, aperture_photometry, CircularAnnulus, RectangularAperture, RectangularAnnulus
from photutils import Background2D, MedianBackground, SExtractorBackground
from photutils.segmentation import detect_sources, SourceFinder, SourceCatalog
### photutils background esitmator
sigma_clip = SigmaClip(sigma = 3.)
bkg_estimator = MedianBackground()
''' Numpy and Scipy modules '''
import numpy as np
from scipy import interpolate, integrate, optimize, stats, ndimage, signal
'''Matplotlib modules'''
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.patches as patches
from matplotlib.lines import Line2D
import matplotlib.animation as animation
%matplotlib inline
'''Some of the following modules are not used, remove if you cannot import'''
## finds all the pathnames matching a specified pattern
import glob
## OS commands
import os
## subprocess with shell commands
import subprocess
## get time stamp
import time
## regular expression
import re
## download packadges? I guess I didn't use these:
# import urllib
# import requests
# import pkg_resources
## Belows are used by my legacy codes on grism simulation, maybe no need any more
# import h5py
# import yaml
# import copy
'''You might need some data from mirage simulator:
https://github.com/spacetelescope/mirage/
'''
# os.environ["MIRAGE_DATA"] = os.environ["MIRAGE_DATA"] + "/mirage_data/"
'''Import CRDS Server to environmental variable'''
os.environ["CDRS_SERVER_URL"]="https://jwst-cdrs.stsci.edu"
## Multiple processing
from multiprocessing import Pool, Lock
## https://github.com/npirzkal/GRISMCONF
import grismconf
## https://github.com/spacetelescope/pysiaf
import pysiaf
## https://pypi.org/project/asdf/
import asdf
## if there is any animation to be display with Jupyter notebook:
from IPython.display import HTML
'''Filter some warning'''
import warnings
warnings.filterwarnings("ignore", message="'obsfix' made the change")
warnings.filterwarnings("ignore", message="'datfix' made the change")
warnings.filterwarnings("ignore", message = "Card is too long")
'''Default nircam zeropoint with native pixel size in LW (0.0629 arcsec/pixel)'''
nircam_LW_ZP = -2.5 * np.log10((u.MJy / u.sr * (0.0629*u.arcsec)**2 / (3631 * u.Jy)).cgs.value)
In [25]:
'''
Matplotlib Drawing:
'''
plt.rcParams.update({'font.size': 16}) ## larger font size
def get_corner_pos(ax, loc = 1, edge = 0.02):
'''get (x, y) for the corner of an axes at given quadrant (1, 2, 3, 4)'''
xlim = ax.get_xlim()
ylim = ax.get_ylim()
xscale, yscale = ax.get_xscale(), ax.get_yscale()
if loc not in [1, 2, 3, 4]: raise ValueError('The parameter `loc` must be 1 or 2 or 3 or 4.')
vec_edge = np.array([1-edge, edge])
vec_x = np.array([[xlim[1], xlim[0], xlim[0], xlim[1]][loc-1], [xlim[0], xlim[1], xlim[1], xlim[0]][loc-1]])
vec_y = np.array([[ylim[1], ylim[1], ylim[0], ylim[0]][loc-1], [ylim[0], ylim[0], ylim[1], ylim[1]][loc-1]])
if xscale == 'linear': x_text = np.dot(vec_edge, vec_x)
else: x_text = 10 ** np.dot(vec_edge, np.log10(vec_x))
if yscale == 'linear': y_text = np.dot(vec_edge, vec_y)
else: y_text = 10 ** np.dot(vec_edge, np.log10(vec_y))
return (x_text, y_text)
def corner_text(ax, s ='', loc = 1, edge = 0.02, **kwargs):
'''Add text at the corner of an axes'''
if loc not in [1, 2, 3, 4]: raise ValueError('The parameter `loc` must be 1 or 2 or 3 or 4.')
ha = ['right', 'left', 'left', 'right']
va = ['top', 'top', 'bottom', 'bottom']
return ax.text(*get_corner_pos(ax, loc, edge), s=s, ha = ha[loc-1], va = va[loc-1], **kwargs)
def rotate(axis, angle):
'''Fundamental rotation matrices.
Rotate by angle measured in degrees, about axis 1 2 or 3'''
if axis not in list(range(1, 4)):
print ('Axis must be in range 1 to 3')
return
r = np.zeros((3, 3))
ax0 = axis-1 #Allow for zero offset numbering
theta = angle * np.pi / 180. # radians(angle)
r[ax0, ax0] = 1.0
ax1 = (ax0+1) % 3
ax2 = (ax0+2) % 3
r[ax1, ax1] = np.cos(theta)
r[ax2, ax2] = np.cos(theta)
r[ax1, ax2] = - np.sin(theta)
r[ax2, ax1] = np.sin(theta)
return r
def attitude(v2, v3, ra, dec, pa):
'''This will make a rotation matrix which rotates a unit vector representing a v2, v3 position
to a unit vector representing an RA, Dec pointing with an assigned position angle
Described in JWST-STScI-001550, SM-12, section 6.1'''
# v2, v3 in arcsec, ra, dec and position angle in degrees
v2d = v2/3600.0
v3d = v3/3600.0
# Get separate rotation matrices
mv2 = rotate(3, -v2d)
mv3 = rotate(2, v3d)
mra = rotate(3, ra)
mdec = rotate(2, -dec)
mpa = rotate(1, -pa)
# Combine as mra*mdec*mpa*mv3*mv2
m = np.dot(mv3, mv2)
m = np.dot(mpa, m)
m = np.dot(mdec, m)
m = np.dot(mra, m)
return m
'''This is the Grism Dispersion function, will be detailed later'''
def fit_disp_order32(data,
a01, a02, a03, a04, a05, a06,
b01, b02, b03, b04, b05, b06,
c01, c02, c03, #c04, c05, c06,
d01, #d02, d03, d04, d05, d06,
):
## data is an numpy array of the shape (3, N)
## - data[0]: x_pixel --> fit with second-degree polynomial
## - data[1]: y_pixel --> fit with second-degree polynomial
## - data[2]: wavelength --> fit with third-degree polynomial
xpix, ypix, dx = data[0] - 1024, data[1] - 1024, data[2] - 3.95
## return dx = dx(x_pixel, y_pixel, lambda)
return ((a01 + (a02 * xpix + a03 * ypix) + (a04 * xpix**2 + a05 * xpix * ypix + a06 * ypix**2)
) +
(b01 + (b02 * xpix + b03 * ypix) + (b04 * xpix**2 + b05 * xpix * ypix + b06 * ypix**2)
) * dx +
(c01 + (c02 * xpix + c03 * ypix) #+ (c04 * xpix**2 + c05 * xpix * ypix + c06 * ypix**2)
) * dx**2 +
(d01 #+ (d02 * xpix + d03 * ypix) + (d04 * xpix**2 + d05 * xpix * ypix + d06 * ypix**2)
) * dx**3
)
func_fit_wave = fit_disp_order32
'''This is the Spectral Tracing function, will be detailed later'''
def fit_disp_order23(data,
a01, a02, a03, a04, a05, a06, a07, a08, a09, a10,
b01, b02, b03, b04, b05, b06, b07, b08, b09, b10,
c01, c02, c03, c04, c05, c06, c07, c08, c09, c10,
):
## data is an numpy array of the shape (3, N)
## - data[0]: x_pixel --> fit with second-degree polynomial
## - data[1]: y_pixel --> fit with second-degree polynomial
## - data[2]: dx --> fit with second-degree polynomial
xpix, ypix, dx = data[0] - 1024, data[1] - 1024, data[2]
## return dy = dy(x_pixel, y_pixel, d_x)
return ((a01 + (a02 * xpix + a03 * ypix) + (a04 * xpix**2 + a05 * xpix * ypix + a06 * ypix**2)
+ (a07 * xpix**3 + a08 * xpix**2 * ypix + a09 * xpix * ypix**2 + a10 * ypix**3)
) +
(b01 + (b02 * xpix + b03 * ypix) + (b04 * xpix**2 + b05 * xpix * ypix + b06 * ypix**2)
+ (b07 * xpix**3 + b08 * xpix**2 * ypix + b09 * xpix * ypix**2 + b10 * ypix**3)
) * dx +
(c01 + (c02 * xpix + c03 * ypix) + (c04 * xpix**2 + c05 * xpix * ypix + c06 * ypix**2)
+ (c07 * xpix**3 + c08 * xpix**2 * ypix + c09 * xpix * ypix**2 + c10 * ypix**3)
) * dx**2
)
def grism_conf_preparation(x0 = 1024, y0 = 1024, pupil = 'R',
fit_opt_fit = np.zeros(30), w_opt = np.zeros(16)):
'''
Prepare grism configuration, dxs, dys, wavelengths based on input (x0, y0) pixel postion
and filter/pupil/module information.
-----------------------------------------------
Parameters
----------
x0, y0 : float
Reference position (i.e., in direct image)
pupil: 'R' or 'C'
pupil of grism ('R' or 'C')
fit_opt_fit: numpy.ndarray, shape: (30,)
polynomial parameters in the perpendicular direction, used by function `fit_disp_order23`
w_opt: numpy.ndarray, shape: (16,)
polynomial parameters in the dispersed direction, used by function `fit_disp_order32`
Returns
-------
dxs, dys : `~numpy.ndarray`
offset of spectral pixel from the direct imaging position
wavs: `~numpy.ndarray`
Array of wavelengths corresponding to dxs and dys
'''
# Load the Grism Configuration file
# GConf = grismconf.Config(os.environ['MIRAGE_DATA'] + "/nircam/GRISM_NIRCAM/V3/" +
# "NIRCAM_%s_mod%s_%s.conf" % (filter, module, pupil))
wave_space = np.arange(2.39, 5.15, 0.01)
disp_space = func_fit_wave(np.vstack((x0 * np.ones_like(wave_space), y0 * np.ones_like(wave_space), wave_space)), *w_opt)
inverse_wave_disp = interpolate.UnivariateSpline(disp_space[np.argsort(disp_space)], wave_space[np.argsort(disp_space)], s = 0, k = 1)
if pupil == 'R':
dxs = np.arange(int(np.min(disp_space)), int(np.max(disp_space))) - x0%1
# Compute wavelength of each of the pixels
wavs = inverse_wave_disp(dxs)
# Compute the dys values for the same pixels
dys = fit_disp_order23(np.vstack((x0* np.ones_like(dxs), y0 * np.ones_like(dxs), dxs)), *fit_opt_fit)
# dxs = np.arange(-1800, 1800, 1) - x0%1
## Compute the t values corresponding to the exact offsets (with old grism config)
# ts = GConf.INVDISPX(order = '+1', x0 = x0, y0 = y0, dx = dxs)
# dys = GConf.DISPY('+1', x0, y0, ts)
# wavs = GConf.DISPL('+1', x0, y0, ts)
# tmp_aper = np.max([0.2, 1.5 * tmp_re_maj * np.cos(tmp_pa), 1.5 * tmp_re_min * np.sin(tmp_pa)])
elif pupil == 'C':
# Compute the dys values for the same pixels
dys = np.arange(int(np.min(disp_space)), int(np.max(disp_space))) - y0%1
# Compute wavelength of each of the pixels
wavs = inverse_wave_disp(dys)
dxs = fit_disp_order23(np.vstack((x0* np.ones_like(dys), y0 * np.ones_like(dys), dys)), *fit_opt_fit)
# dys = np.arange(-1800, 1800, 1) - y0%1
## Compute the t values corresponding to the exact offsets (with old grism config)
# ts = GConf.INVDISPY(order = '+1', x0 = x0, y0 = y0, dy = dys)
# dxs = GConf.DISPX('+1', x0, y0, ts)
# wavs = GConf.DISPL('+1', x0, y0, ts)
return (dxs, dys, wavs)
'''Lienar Function'''
linear = lambda x, k, b: x * k + b
'''gaussian function'''
gauss = lambda x, x0, flux, FWHM : (flux / FWHM / 1.064467) * np.exp(-(x - x0)**2 / (2 * (FWHM/2.354820)**2))
'''gaussian function + continuum'''
gauss_cont_prof = lambda x, x0, flux, FWHM, k, b : (flux / FWHM / 1.064467) * np.exp(-(x - x0)**2 / (2 * (FWHM/2.354820)**2)) + k * (x - 1024) + b
Not all of the source in your input catalog could yield spectra on your 2D grism image. For example, below is the spectral coverage for F356W grism-R in module A and B:
In this plot, everything is in the detector frame (X and Y; e.g., the detector range is X=0-2047, Y=0-2047 which is indicated by the red box). Based on pre-launch estimate, the regions that can be perfectly picked off by the POM are shown in shallow grey color, and the coronagraph regions are in white - usally the spectral quality will be affected if your source falls in these regions. Only sources in the blue (module-A) and orange (module-B) can yield complete spectra with the current filter setting (F356W). The pre-launch estimate of full spectral coverage region is indicated by the black boxes, so you can see some differences. The yielded spectra alongside with the dispersion directions are indicated by the arrows. If the source falls in the hatched region, it could yield partial spectral coverage. You can get access to all of the spectral coverage plots of all NIRCam LW filters and grisms at http://u.arizona.edu/~fengwusun/data/spectral_cov/.
The following codes are used to check whether the source could enter Pick-Off Mirror (POM) and yield spectrum:¶
In [26]:
## module A/B LW POM transmission:
## These files can be downloaded at: https://github.com/npirzkal/GRISM_NIRCAM/tree/master/V4
modA_LW_POM_fits = fits.open(os.environ['MIRAGE_DATA'] + '/nircam/GRISM_NIRCAM/current/NIRCAM_LW_POM_ModA_trans.fits')
modB_LW_POM_fits = fits.open(os.environ['MIRAGE_DATA'] + '/nircam/GRISM_NIRCAM/current/NIRCAM_LW_POM_ModB_trans.fits')
modA_LW_POM_trans = modA_LW_POM_fits[1].data # np.clip(np.int8(modA_LW_POM_fits[1].data), 0, 1)
modB_LW_POM_trans = modB_LW_POM_fits[1].data # np.clip(np.int8(modB_LW_POM_fits[1].data), 0, 1)
xy_start_POM_modA = modA_LW_POM_fits[1].header['NOMXSTRT'], modA_LW_POM_fits[1].header['NOMYSTRT']
xy_start_POM_modB = modB_LW_POM_fits[1].header['NOMXSTRT'], modA_LW_POM_fits[1].header['NOMYSTRT']
def is_pickoff(x, y, module = 'A'):
'''
Input x/y pixel(s) and module, find whether the source will be picked off by NIRCam POM.
'''
### get module and start point
if module.upper() == 'A': xy_start, POM_trans = xy_start_POM_modA, modA_LW_POM_trans
elif module.upper() == 'B': xy_start, POM_trans = xy_start_POM_modB, modB_LW_POM_trans
else: raise ValueError("Only \"A\" or \"B\" is allowed for module.")
if type(x) == np.ndarray: tmp_x, tmp_y = np.int32(x), np.int32(y)
elif type(x) == list: tmp_x, tmp_y = np.array(x, dtype = np.int32), np.array(y, dtype = np.int32)
elif type(x) == int: tmp_x, tmp_y = np.array([x], dtype = np.int32), np.array([y], dtype = np.int32)
elif type(x) == float: tmp_x, tmp_y = np.array([x], dtype = np.int32), np.array([y], dtype = np.int32)
else: raise TypeError("Only integers/floats or their list/np.array() are allowed.")
if len(tmp_x) != len(tmp_y): raise ValueError("x, y list should be of the same length.")
tmp_x = tmp_x + xy_start[0]
tmp_y = tmp_y + xy_start[1]
arr_trans = []
for j in range(len(tmp_x)):
if (tmp_x[j] < 0)|(tmp_y[j]<0)|(tmp_x[j] > POM_trans.shape[1]-1)|(tmp_y[j] > POM_trans.shape[0]-1):
arr_trans.append(0)
else:
arr_trans.append(POM_trans[tmp_x[j], tmp_x[j]])
return np.array(arr_trans, dtype = int)
def is_pickoff_PS(x, y, module = 'A', filter = 'F322W2', pupil = 'C'):
'''
Input x/y pixel(s) and module,
find whether the source will be picked off by NIRCam POM and yield partial/complete spectra.
'''
### Get module and start point
## partial spectral coverage region fits files (`path_PS1) can be downloaded at:
## http://gxn.as.arizona.edu/~sunfengwu/spec_cov/
## now available for both wide and medium LW filter now.
path_PS = '/groups/egami/mirage_data/mirage_data/nircam/GRISM_NIRCAM/FSun_SpecCov_%s_%s_%s.fits' % (filter, module, pupil)
path_POM = os.environ['MIRAGE_DATA'] + '/nircam/GRISM_NIRCAM/current/NIRCAM_LW_POM_Mod%s_trans.fits' % module
if os.path.isfile(path_PS) == False:
raise FileNotFoundError("Spectral Coverage file %s not found!" % path_PS)
# if os.path.isfile(path_POM) == False:
# raise FileNotFoundError("POM Coverage file %s not found!" %path_POM)
spec_cov_fits = fits.open(path_PS)
xy_start = (spec_cov_fits[0].header['NOMXSTRT'], spec_cov_fits[0].header['NOMYSTRT'])
POM_trans = spec_cov_fits[0].data # * np.int16(fits.getdata(path_POM) == 1)
spec_cov_fits.close()
if type(x) == np.ndarray: tmp_x, tmp_y = np.int32(x), np.int32(y)
elif type(x) == list: tmp_x, tmp_y = np.array(x, dtype = np.int32), np.array(y, dtype = np.int32)
elif type(x) == int: tmp_x, tmp_y = np.array([x], dtype = np.int32), np.array([y], dtype = np.int32)
elif type(x) == float: tmp_x, tmp_y = np.array([x], dtype = np.int32), np.array([y], dtype = np.int32)
else: raise TypeError("Only integers/floats or their list/np.array() are allowed.")
if len(tmp_x) != len(tmp_y): raise ValueError("x, y list should be of the same length.")
tmp_x = tmp_x + xy_start[0]
tmp_y = tmp_y + xy_start[1]
arr_trans = []
for j in range(len(tmp_x)):
if (tmp_x[j] < 0)|(tmp_y[j]<0)|(tmp_x[j] > POM_trans.shape[1]-1)|(tmp_y[j] > POM_trans.shape[0]-1):
arr_trans.append(0)
### For now, I don't trust the spectra from Coronagraphy region (y >~ 2250);
### I'm masking all of these frames:
elif tmp_y[j] > 2250:
arr_trans.append(0)
else:
arr_trans.append(POM_trans[tmp_y[j], tmp_x[j]])
return np.array(arr_trans, dtype = float)
2. Get list of RATE files (stage-1 products; from MAST)¶
After downloading your rate files from MAST https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html or reducing *_uncal.fits to *_rate.fits files, you can get a list of the rate files and start analyzing the observational setups (e.g., filter, module, pupil, targets) of these observations:
In [13]:
'''
List of the rate files on your local directory:
Here I'm only selecting module A/B long wavelength images by wildcards: `jw*nrc[a-b]long_rate`
'''
# list_rate = np.array(glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/all_rate/obs24_F322W2_C/jw*nrc[a-b]long_rate.fits'))
list_rate = np.array(glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/all_rate/*/jw*nrc[a-b]long_rate.fits'))
list_rate.sort()
'''List of filter, pupil, module, target, exposure time, etc.'''
list_filter = []
list_pupil = []
list_module = []
list_target = []
list_exp = []
for x in list_rate:
### read fits header one by one:
tmp_hd = fits.getheader(x)
list_filter.append(tmp_hd['FILTER'])
list_pupil.append(tmp_hd['PUPIL'])
list_module.append(tmp_hd['MODULE'])
list_target.append(tmp_hd['TARGPROP'])
list_exp.append(tmp_hd['EFFEXPTM'])
list_filter = np.array(list_filter)
list_pupil = np.array(list_pupil)
list_module = np.array(list_module)
list_target = np.array(list_target)
list_exp = np.array(list_exp)
print(np.unique(list_filter), '\n',
np.unique(list_pupil), '\n',
np.unique(list_module), '\n',
np.unique(list_target))
Now at above, you can find out what are the exposures taken in the LW grism (look at the filter and pupil, pupil should be GRISMR or GRISMC). You can also check the number of available exposures using the following command:
In [14]:
### check number of grism exposures
np.sum((list_filter == 'F322W2') & (list_pupil == 'GRISMC') & (list_module == 'A'))
Out[14]:
3. Reduce Grism exposures to stage level-1.5 (i.e., bkg subtracted, wcs assigned)¶
In the following steps, we will reduce the grism data to level 1.5. This includes:
- World Coordinate System (WCS) assginment.
With the rate files you get above, it usually doesn't contain WCS in the sci fits header. You will need to add WCS information partially using the STScI stage-2 pipeline. The WCS information in such images will be use to convert (RA, Dec) coordinates of your sources to (x, y) arrays assuming the pupil wheel is clear.
- Flat Fielding
Currently, there is no available NIRCam grism flat field data. One good approximation is to use the image flat field data constructed in-flight, which is quite accurate (guestimate accuracy ~3% per pixel).
- Background Subtraction
Currently, the theoretical sky background for NIRCam grism observations are all pre-launch estimate and they are not accurate enough. Therefore, I suggest users to construct super sky background using all flat-fielded exposure taken with the same filter+pupil+module combination, and subtract the background there. We might still end up with some background residuals, which could be subtracted using 2D background modeling method.
- 1/f Noise Subtraction
I'm using the code tshirt/ROEBA which can subtract the horizontal 1/f noise stripe for C grism data. However, right now I don't have good idea on how to deal with R grism data where the dispersion direction is the same as the readout direction, so I'm skipping it.
Below is a plot showing the steps of flat-fielding, background subraction and 1/f noise subtraction on IceAge data (F322W2, module B, grism C):
In [27]:
'''
## Stage 2.1 Background subtraction + hot pixel subtraction
## 2.1.1 get list of rate.fits by module and pupil:
'''
tmp_filter = 'F444W' # np.unique(list_filter)[0]
list_rate_this_band = list_rate[(list_filter == tmp_filter) & (list_pupil != 'CLEAR')]
print(len(list_rate_this_band), 'files,\n', list_rate_this_band[0])
In [28]:
'''Python modules to run JWST pipelines & 1/f noise subtractions'''
from jwst.pipeline import Detector1Pipeline
from jwst import assign_wcs, datamodels, tweakreg
from jwst.flatfield import flat_field
## https://github.com/spacetelescope/crds
import crds
## https://tshirt.readthedocs.io/en/latest/index.html
from tshirt.pipeline.instrument_specific import rowamp_sub
'''Directory of grism exposures for output -> Change it as your output folder'''
calibrated_dir = '/xdisk/egami/fengwusun/ERSs/IceAge/grism/'
if os.path.isdir(calibrated_dir + 'plots/') == False:
os.mkdir(calibrated_dir + 'plots/') # used to save analytic plots
def background_grism_stage2(arr_grism_rate, tmp_module, tmp_pupil, tmp_filter):
'''
Function to creat Median background for each module & Pupil
'''
print('\n>> Create Median Background for %s - module %s - grism %s' % (tmp_filter, tmp_module, tmp_pupil))
tmp_med_bkg = sigma_clipped_stats(arr_grism_rate, axis = -1, sigma_upper = 2.0, maxiters = 10)[1]
tmp_med_bkg_path = calibrated_dir + 'median_bkg_%s_mod%s_%s.fits' % (tmp_filter, tmp_module, tmp_pupil)
tmp_med_bkg_fits = fits.HDUList(fits.PrimaryHDU(tmp_med_bkg))
tmp_med_bkg_fits[0].header['object'] = 'BKG_%s_mod%s_%s' % (tmp_filter, tmp_module, tmp_pupil)
tmp_med_bkg_fits[0].header['filter'] = tmp_filter
tmp_med_bkg_fits[0].header['module'] = tmp_module
tmp_med_bkg_fits[0].header['pupil'] = tmp_pupil
tmp_med_bkg_fits.writeto(tmp_med_bkg_path, overwrite = True)
print('>> Saved background models to %s' % tmp_med_bkg_path)
def assignwcs_grism_stage2(rate_grism_file):
'''
Function to reduce grism data to level 1.5 (subtract constructed bkg; assign wcs)
'''
rate_grism_fits = fits.open(rate_grism_file)
tmp_grism_hd = rate_grism_fits[0].header
### read distortion assignment
siaf_file = crds.getreferences(tmp_grism_hd, reftypes = ['distortion'], ignore_cache = False)['distortion']
### find available background
tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil']
print(' read %s' % tmp_path_rate, '(%s - %s - %s)' % (tmp_filter, tmp_module, tmp_pupil))
### correct WCS and background (need to fool STScI pipeline to assume it as NRC image)
grism_wcs_step = assign_wcs.assign_wcs_step.AssignWcsStep(override_distortion = siaf_file)
rate_grism_fits[0].header['EXP_TYPE'] = 'NRC_IMAGE'
grism_image = datamodels.image.ImageModel(rate_grism_fits)
'''don\'t subtract background here'''
# try:
# tmp_bkg_img = fits.getdata(calibrated_dir + 'median_bkg_%s_mod%s_%s.fits' % (tmp_filter, tmp_module, tmp_pupil))
# grism_image.data = np.float32(grism_image.data - tmp_bkg_img)
# except FileNotFoundError:
# pass
grism_with_wcs = grism_wcs_step(grism_image)
### correct for flat:
tmp_flat_path = crds.getreferences({'INSTRUME': tmp_grism_hd['INSTRUME'], 'READPATT': tmp_grism_hd['READPATT'],
'SUBARRAY': tmp_grism_hd['SUBARRAY'],
'DATE-OBS': tmp_grism_hd['DATE'].split('T')[0],'TIME-OBS': tmp_grism_hd['DATE'].split('T')[1],
'DETECTOR': tmp_grism_hd['DETECTOR'], 'CHANNEL': tmp_grism_hd['CHANNEL'],
'MODULE': tmp_grism_hd['MODULE'], 'EXP_TYPE': 'NRC_IMAGE', #'EXP_TYPE': 'NRC_WFSS',
'FILTER': tmp_grism_hd['FILTER'], 'PUPIL': 'CLEAR', # 'PUPIL': 'GRISMR', # =
})['flat']
### error of Flat
# 1 / np.percentile(fits.getdata(tmp_flat_path, 'err')/ fits.getdata(tmp_flat_path, 'sci'), [16, 50, 84])
flat_field.do_flat_field(grism_with_wcs, datamodels.FlatModel(tmp_flat_path))
grism_save_path = rate_grism_file.replace('rate.fits', 'rate_lv1.5.fits')
grism_save_path = calibrated_dir + grism_save_path.split('/')[-1]
print("Saving file with proper WCS in", grism_save_path)
grism_with_wcs.save(grism_save_path)
def get_crds_dict_from_fits_header(tmp_hd):
'''get the CRDS directory from the fits header'''
crds_dict = {}
crds_dict['INSTRUME'] = tmp_hd['INSTRUME'].upper()
crds_dict['READPATT'], crds_dict['SUBARRAY'] = tmp_hd['readpatt'].upper(), tmp_hd['SUBARRAY'].upper()
crds_dict['DATE-OBS'], crds_dict['TIME-OBS'] = tmp_hd['DATE-OBS'], tmp_hd['TIME-OBS']
crds_dict['DETECTOR'] = tmp_hd['DETECTOR']
if crds_dict['INSTRUME'] == 'NIRCAM':
if crds_dict['DETECTOR'] in ['NRCALONG', 'NRCBLONG']: crds_dict['CHANNEL'] = 'LONG'
else: crds_dict['CHANNEL'] = 'SHORT'
crds_dict['MODULE'] = tmp_hd['module']
crds_dict['EXP_TYPE'], crds_dict['FILTER'], crds_dict['PUPIL'] = tmp_hd['EXP_TYPE'], tmp_hd['filter'], tmp_hd['pupil']
return crds_dict
In [29]:
'''
## Stage 2.2.a WCS corrections + Flat Fielding subtraction
'''
n_procs = 8
### Uncomment the following lines to do this step in parallel
# with Pool(np.min([n_procs, len(list_rate_this_band)])) as pool:
# pool.map(assignwcs_grism_stage2, list_rate_this_band[:])
### after this step, change the suffix of the files to `rate_lv1.5.fits`
all_rate_list = np.array([calibrated_dir + os.path.basename(x).replace('rate.fits', 'rate_lv1.5.fits')
for x in list_rate_this_band])
In [ ]:
'''Check image of flat-fielded grism exposure'''
a_c_sci = fits.getdata(all_rate_list[0])
plt.close()
plt.figure(figsize = (10, 10))
plt.imshow(a_c_sci,
vmin = np.nanpercentile(a_c_sci, 0.5), vmax = np.nanpercentile(a_c_sci, 99.5), origin = 'lower')
plt.tight_layout()
In [763]:
'''
## Stage 2.2.b Generate Super Sky Background
'''
print("- Background + Hot pixel construction for %s - " % tmp_filter)
count_ac, count_bc, count_ar, count_br = 0, 0, 0, 0
for i, tmp_path_rate in enumerate(all_rate_list):
tmp_grism_hd = fits.getheader(tmp_path_rate)
tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil']
print('>>> [%3d] ' % i, 'read %s' % tmp_path_rate,
'(%s %s %s)' % (tmp_filter, tmp_module, tmp_pupil))
tmp_grism_rate = fits.getdata(tmp_path_rate, 'sci')
if (tmp_module == 'A') & (tmp_pupil == 'GRISMC'):
if count_ac == 0: arr_grism_rate_ac = tmp_grism_rate
else: arr_grism_rate_ac = np.dstack((arr_grism_rate_ac, tmp_grism_rate))
count_ac += 1
if (tmp_module == 'B') & (tmp_pupil == 'GRISMC'):
if count_bc == 0: arr_grism_rate_bc = tmp_grism_rate
else: arr_grism_rate_bc = np.dstack((arr_grism_rate_bc, tmp_grism_rate))
count_bc += 1
if (tmp_module == 'A') & (tmp_pupil == 'GRISMR'):
if count_ar == 0: arr_grism_rate_ar = tmp_grism_rate
else: arr_grism_rate_ar = np.dstack((arr_grism_rate_ar, tmp_grism_rate))
count_ar += 1
if (tmp_module == 'B') & (tmp_pupil == 'GRISMR'):
if count_br == 0: arr_grism_rate_br = tmp_grism_rate
else: arr_grism_rate_br = np.dstack((arr_grism_rate_br, tmp_grism_rate))
count_br += 1
print(count_ar, count_ac, count_br, count_bc)
'''
## generating background by module and pupil:
'''
### if you only have C grism data, use below:
# list_arr_grism_rate = [arr_grism_rate_ac, arr_grism_rate_bc]
# list_arr_grism_module = ['A', 'B']
# list_arr_grism_pupil = ['GRISMC', 'GRISMC']
# list_arr_grism_filter = [tmp_filter] * 2
### if you only have R grism data, use below:
# list_arr_grism_rate = [arr_grism_rate_ar, arr_grism_rate_br]
# list_arr_grism_module = ['A', 'B']
# list_arr_grism_pupil = ['GRISMR', 'GRISMR']
# list_arr_grism_filter = [tmp_filter] * 2
### if you have both R and C grism data, use below:
list_arr_grism_rate = [arr_grism_rate_ar, arr_grism_rate_ac, arr_grism_rate_br, arr_grism_rate_bc]
list_arr_grism_module = ['A', 'A', 'B', 'B']
list_arr_grism_pupil = ['GRISMR', 'GRISMC', 'GRISMR', 'GRISMC']
list_arr_grism_filter = [tmp_filter] * 4
n_procs = 8
### Uncomment the following lines to do this step in parallel
# with Pool(np.min([n_procs, len(all_rate_list)])) as pool:
# pool.starmap(background_grism_stage2, zip(list_arr_grism_rate, list_arr_grism_module,
# list_arr_grism_pupil, list_arr_grism_filter))
In [ ]:
'''Check the grism image:'''
tmp_i = np.random.randint(len(all_rate_list))
tmp_i = 51
print(all_rate_list[tmp_i])
tmp_grism_img = fits.getdata(all_rate_list[tmp_i])
tmp_grism_hd = fits.getheader(all_rate_list[tmp_i])
tmp_grism_bkg = fits.getdata(calibrated_dir + 'median_bkg_%s_mod%s_%s.fits' %
(tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil']))
fig, ax = plt.subplots(2, 2, figsize = (15, 15))
ax = ax.flatten()
'''ax-0: flat-fielded grism image'''
ax[0].imshow(tmp_grism_img, origin = 'lower',
vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[0], s = '(1) Flat-Fielded', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
'''ax-1: try to subtract sigma-clipped median sky background'''
if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'):
pass ## this part of IceAge data is affected by a few very bright stars
else: tmp_grism_img = tmp_grism_img - tmp_grism_bkg
ax[1].imshow(tmp_grism_img, origin = 'lower',
vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[1], s = '(2) Median Sky BKG Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'):
corner_text(ax[1], s = '(Not applied for modB - GRISMR)', color = 'w', loc = 4, weight = 'bold', fontsize = 20)
# mask-1 for background / 1-over-f noise subtraction: from segment map
_, tmp_med, tmp_rms = sigma_clipped_stats(tmp_grism_img[100:700,100:700])
segment_map = detect_sources(tmp_grism_img - tmp_med, tmp_rms * 2.0, npixels = 100)
### further background subtraction?
sigma_clip = SigmaClip(sigma = 2.5)
bkg_estimator = SExtractorBackground()
if (tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('B', 'R'):
tmp_mask = segment_map.data!=0
else: tmp_mask = np.zeros_like(segment_map.data)
bkg = Background2D(tmp_grism_img, (24, 24), filter_size = (5, 5), mask = tmp_mask,
sigma_clip = sigma_clip, bkg_estimator = bkg_estimator)
tmp_grism_img = tmp_grism_img - bkg.background
'''ax-2: try to subtract modeled 2D sky background'''
ax[2].imshow(tmp_grism_img, origin = 'lower',
vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[2], s = '(3) Modeled 2D BKG Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
### for grism C: subtract 1/f noise
if tmp_grism_hd['pupil'][-1] == 'C':
_, tmp_med, tmp_rms = sigma_clipped_stats(tmp_grism_img[100:1700,100:1700].flatten()[::10])
segment_map = detect_sources(tmp_grism_img, tmp_rms * 5.0, npixels = 20)
rowSub, modelImg = rowamp_sub.do_backsub(tmp_grism_img,
amplifiers = 4, backgMask = segment_map.data == 0)
tmp_grism_img = rowSub
'''ax-3: try to subtract 1/f noise'''
ax[3].imshow(tmp_grism_img, origin = 'lower',
vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[3], s = '(4) 1/f Noise Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
if tmp_grism_hd['pupil'][-1] == 'R':
corner_text(ax[3], s = '(Not applied for GRISMR)', color = 'w', loc = 4, weight = 'bold', fontsize = 20)
tmp_title = '%s (%s - %s - %s)' % (os.path.basename(all_rate_list[tmp_i]).replace('_lv1.5.fits', ''),
tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'])
fig.text(0.5, 0.985, s = tmp_title, ha = 'center', va = 'top', fontsize = 24,)
plt.subplots_adjust(left = 0.02, right = 0.98, top = 0.95, bottom = 0.01, wspace = 0.06, hspace = 0.04)
for tmp_ax in ax:
tmp_ax.set_xticks([]); tmp_ax.set_yticks([])
# plt.savefig(calibrated_dir + 'plots/%s_mod%s_%s_%2d_%s.pdf' % (tmp_grism_hd['filter'], tmp_grism_hd['module'],
# tmp_grism_hd['pupil'], tmp_i,
# tmp_title.split(' ')[0]), dpi = 100)
plt.show()
In [ ]:
'''If necessary, check background / segmentation / 1-over-f noise model'''
plt.figure(figsize = (8, 8))
plt.imshow(bkg.background, origin = 'lower', vmax = np.nanpercentile(bkg.background, 99.5))
print(np.nanpercentile(bkg.background, [0, 100]))
plt.imshow(tmp_mask, origin = 'lower')
# plt.imshow(modelImg)
In [22]:
'''
## Stage 2.2.c Background Ground (1/f noise subtraction) - For loop or pool multi-processing
'''
def my_grism_bkg_subtraction(rate_grism_file):
print('>>> subtract background for %s ' % rate_grism_file)
tmp_grism_fits = fits.open(rate_grism_file)
tmp_grism_img = tmp_grism_fits[1].data
tmp_grism_hd = tmp_grism_fits[0].header
tmp_grism_bkg = fits.getdata(calibrated_dir + 'median_bkg_%s_mod%s_%s.fits' %
(tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil']))
plt.close()
fig, ax = plt.subplots(2, 2, figsize = (15, 15))
ax = ax.flatten()
'''ax-0: flat-fielded grism image'''
ax[0].imshow(tmp_grism_img, origin = 'lower',
vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[0], s = '(1) Flat-Fielded', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
'''ax-1: subtract sigma-clipped median sky background'''
if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'): pass
else:
tmp_grism_img = tmp_grism_img - tmp_grism_bkg
tmp_grism_fits[0].header['HISTORY'] = 'sigma-clipped median sky background subtracted by F. Sun on %s' % time.strftime("%Y/%m/%d", time.localtime())
ax[1].imshow(tmp_grism_img, origin = 'lower',
vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[1], s = '(2) Median Sky BKG Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'):
corner_text(ax[1], s = '(Not applied for modB - GRISMR)', color = 'w', loc = 4, weight = 'bold', fontsize = 20)
# mask-1: from segment map
_, tmp_med, tmp_rms = sigma_clipped_stats(tmp_grism_img[100:700,100:700])
segment_map = detect_sources(tmp_grism_img - tmp_med, tmp_rms * 1.0, npixels = 100)
### further background subtraction?
sigma_clip = SigmaClip(sigma = 2.5)
bkg_estimator = SExtractorBackground()
if (tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]) == ('F322W2', 'B', 'R'):
tmp_mask = segment_map.data!=0
# bkg = Background2D(tmp_grism_img, (36, 36), filter_size = (5, 5), mask = tmp_mask,
# sigma_clip = sigma_clip, bkg_estimator = bkg_estimator, exclude_percentile=50.)
else:
tmp_mask = np.zeros_like(segment_map.data)
bkg = Background2D(tmp_grism_img, (24, 24), filter_size = (5, 5), mask = tmp_mask,
sigma_clip = sigma_clip, bkg_estimator = bkg_estimator, exclude_percentile=50.)
tmp_grism_img = tmp_grism_img - bkg.background
tmp_grism_fits[0].header['HISTORY'] = 'SExtractor-modeled 2D background subtracted by F. Sun on %s' % time.strftime("%Y/%m/%d", time.localtime())
'''ax-2: subtract modeled 2D sky background'''
ax[2].imshow(tmp_grism_img, origin = 'lower',
vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[2], s = '(3) Modeled 2D BKG Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
### for grism C: subtract 1/f noise
if tmp_grism_hd['pupil'][-1] == 'C':
_, tmp_med, tmp_rms = sigma_clipped_stats(tmp_grism_img[100:1700,100:1700].flatten()[::10])
segment_map = detect_sources(tmp_grism_img, tmp_rms * 5.0, npixels = 20)
rowSub, modelImg = rowamp_sub.do_backsub(tmp_grism_img,
amplifiers = 4, backgMask = segment_map.data == 0)
tmp_grism_img = rowSub
tmp_grism_fits[0].header['HISTORY'] = '1/f noise subtracted by F. Sun on %s' % time.strftime("%Y/%m/%d", time.localtime())
'''ax-3: subtract 1/f noise'''
ax[3].imshow(tmp_grism_img, origin = 'lower',
vmin = np.nanpercentile(tmp_grism_img, 2.5), vmax = np.nanpercentile(tmp_grism_img, 97.5))
corner_text(ax[3], s = '(4) 1/f Noise Subtracted', color = 'w', loc = 2, weight = 'bold', fontsize = 20)
if tmp_grism_hd['pupil'][-1] == 'R':
corner_text(ax[3], s = '(Not applied for GRISMR)', color = 'w', loc = 4, weight = 'bold', fontsize = 20)
tmp_title = '%s (%s - %s - %s)' % (os.path.basename(rate_grism_file).replace('_lv1.5.fits', ''),
tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'])
fig.text(0.5, 0.985, s = tmp_title, ha = 'center', va = 'top', fontsize = 24,)
plt.subplots_adjust(left = 0.02, right = 0.98, top = 0.95, bottom = 0.01, wspace = 0.06, hspace = 0.04)
for tmp_ax in ax:
tmp_ax.set_xticks([]); tmp_ax.set_yticks([])
plt.savefig(calibrated_dir + 'plots/%s_mod%s_%s_%s.pdf' % (tmp_grism_hd['filter'], tmp_grism_hd['module'],
tmp_grism_hd['pupil'], tmp_title.split(' ')[0]), dpi = 100)
plt.close()
'''Write Fits File'''
tmp_grism_fits[1].data = tmp_grism_img
tmp_grism_fits.writeto(rate_grism_file, overwrite = True)
tmp_grism_fits.close()
n_procs = 8
### Uncomment the following lines to do this step in parallel:
# with Pool(np.min([n_procs, len(all_rate_list)])) as pool:
# pool.map(my_grism_bkg_subtraction, all_rate_list)
### Or, uncomment the following lines to do this step in loop:
# for i, x in enumerate(all_rate_list):
# print(i); my_grism_bkg_subtraction(x)
At this step, you should have completed the stage-2a reduction for the LW grism images. The files (jw*nrc[a-b]long_rate_lv1.5.fits) you get are flat-fielded, wcs-assigned and background subtracted LW grism images, almost ready for spectral extraction. However, there is one more thing that you will need to pay attention, which is astrometry.
The spectral and spatial resolution of NIRCam (grism) is so high that a small astrometric offset (0.3" -> 5 LW pix) can give you a completely trash 2D or 1D spectra. The quality of astrometric correction usually limits the quality of spectral extraction and flux/wavelength calibration. One good news is that we have SW simultaneous direct imaging data which can be used to empirically calibrate the LW grism astrometry. This is, assuming that when SW direct imaging data is registered to certain astrometric catalog by applying dRA, dDEc and theta (i.e., x/y offset and rotation), the same offset will be applied to the LW grism fits header and the LW grism astrometry should be good. Such a registration relies on the internal alignment between NIRCam SW and LW instrument aperture and it may introduce astrometric residuals. However, these residuals should be stable and have been included in the grism spectral tracing models.
When using the different SW/LW filter combinations, different filters will also introduce astrometric offset called "filter offset". In my personal opinion, this has not been well understood and calibrated throughout commissioning and it can cause most of the tracing and wavelength error (a few pixels). We will get to this point later.
4. Astrometry Characterization using SW Direct imaging¶
In [32]:
'''Getting all SW rate fits files by wildcards'''
list_rate_sw = np.array(glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/all_rate/*/jw*_nrc[a-b][1-4]_rate.fits'))
list_rate_sw.sort()
print(len(list_rate_sw))
'''List of filter, pupil, module, target, exposure time, etc.'''
list_sw_filter = []
list_sw_pupil = []
list_sw_module = []
list_sw_target = []
list_sw_exp = []
for x in list_rate_sw:
tmp_hd = fits.getheader(x)
list_sw_filter.append(tmp_hd['FILTER'])
list_sw_pupil.append(tmp_hd['PUPIL'])
list_sw_module.append(tmp_hd['MODULE'])
list_sw_target.append(tmp_hd['TARGPROP'])
list_sw_exp.append(tmp_hd['EFFEXPTM'])
list_sw_filter = np.array(list_sw_filter)
list_sw_pupil = np.array(list_sw_pupil)
list_sw_module = np.array(list_sw_module)
list_sw_target = np.array(list_sw_target)
### print SW filters and modules
print(np.unique(list_sw_filter), np.unique(list_sw_module), np.unique(list_sw_target))
For IceAge data in our example, if I'm looking at the F444W grism data, the SW simultaneous exposures are always taken in F200W filter, so I'm only selecting SW rate files with list_sw_filter = 'F200W' . For these data, I have run the following steps in another jupyter notebook:
reduced all
*_rate.fitsfiles to stage-2 as*_cal_fits, including flat-fielding, distortion correction, photometric zeropoint assignment and 1/f noise subtraction.register the WCS of all
*_cal_fitsto an astrometric catalog, constructed at LW.
Therefore, it might be good to detect stars in stage-2 products *_cal_fits, but the WCS information there has been corrected. What we will need is catalogs of stars in all SW images without any astrometric correction, so we can correct for x/y offset and rotration ourselves and pass that to LW. If you haven't corrected anything for _rate.fits and _cal.fits files, great, you can directly run daofind() on the _cal.fits files, get catalogs and compute the astrometric offsets. If you have corrected the astrometry (like me, I did SW/LW imaging first for all IceAge data), we will restore the original astrometry to extract the catalog. Therefore, for the following code block, please modify the scripts to let it suit your cases.
In [34]:
all_rate_sw_list = list_rate_sw[list_sw_filter == 'F200W']
print(len(all_rate_sw_list), 'SW rate files selected.')
'''Path to my direct imaging folder
which has rate files, wcs uncorrected (original),
and cal files, wcs corrected (modified).
'''
direct_image_dir = '/xdisk/egami/fengwusun/ERSs/IceAge/direct/'
# i = 72
# tmp_rate_sw = all_rate_sw_list[i]
def my_daofind_sw_fits(tmp_rate_sw):
'''
Since SW images have been astrometrically corrected and passed through stage-2 pipelines,
Here we only restore the original atrometry and use the images provided by stage-2 pipeline
'''
## assign wcs
rate_sw_fits = fits.open(tmp_rate_sw)
rate_sw_hd = rate_sw_fits[0].header
### read distortion assignment
siaf_file = crds.getreferences(rate_sw_hd, reftypes = ['distortion'], ignore_cache = False)['distortion']
### assign default WCS
rate_sw_wcs_step = assign_wcs.assign_wcs_step.AssignWcsStep(override_distortion = siaf_file)
rate_sw_IM = datamodels.image.ImageModel(rate_sw_fits)
rate_sw_IM_wcs = rate_sw_wcs_step(rate_sw_IM)
tmp_wcs = rate_sw_IM_wcs.get_fits_wcs() ## this is the default wcs from observational parameters
rate_sw_IM_wcs.close()
rate_sw_fits.close()
### well-processed SW cal.fits files
direct_sw_fits = fits.open(direct_image_dir + os.path.basename(tmp_rate_sw).replace('_rate.fits', '_cal.fits'))
tmp_coord_center = wcs.utils.pixel_to_skycoord(1024.5, 1024.5, tmp_wcs)
### construct detection image
tmp_detect_img = np.nan_to_num(direct_sw_fits['sci'].data)
bkg_estimator = MedianBackground()
bkg = Background2D(tmp_detect_img, (64, 64), filter_size = (5, 5),
sigma_clip = SigmaClip(sigma = 3.), bkg_estimator = bkg_estimator)
tmp_detect_std = sigma_clipped_stats(tmp_detect_img[100:1500,100:1500].flatten()[::7],
sigma = 3, maxiters = 10)[-1]
tmp_detect_img = (tmp_detect_img - bkg.background) / tmp_detect_std
# tmp_detect_img[direct_sw_fits['dq'].data >= 1] = 0
### run DAOStarFinder - change parameters (FWHM or Threshold if needed)
daofind = DAOStarFinder(fwhm = 2.5, threshold = 8.)
tmp_tb_daofind = daofind(tmp_detect_img)
tmp_tb_daofind['detector'] = direct_sw_fits[0].header['detector'].lower()
tmp_tb_daofind['skycoord'] = wcs.utils.pixel_to_skycoord(tmp_tb_daofind['xcentroid'], tmp_tb_daofind['ycentroid'], tmp_wcs)
### Save SW catalogs
ascii.write(tmp_tb_daofind, calibrated_dir + 'astrom/%s' % os.path.basename(tmp_rate_sw).replace('_rate.fits', '_daofind.cat'),
format = 'ecsv', overwrite = True)
direct_sw_fits.close()
if os.path.isdir(calibrated_dir + 'astrom/') == False:
os.mkdir(calibrated_dir + 'astrom/') # used to save astrometry-related files
n_procs = 8
### Uncomment the following lines to generate source catalogs in parallel:
# with Pool(np.min([n_procs, len(all_rate_sw_list)])) as pool:
# pool.map(my_daofind_sw_fits, all_rate_sw_list)
'''If some cores are corrupted, run the following for unfinished fits'''
idx_unfinished = np.array([os.path.isfile(calibrated_dir + 'astrom/%s' % os.path.basename(x).replace('_rate.fits', '_daofind.cat')) == False
for x in all_rate_sw_list])
print(np.sum(idx_unfinished), 'unfinished.')
# with Pool(np.min([n_procs, len(all_rate_sw_list[idx_unfinished])])) as pool:
# pool.map(my_daofind_sw_fits, all_rate_sw_list[idx_unfinished])
With the SW catalogs produced above, you can then use an astrometric catalog to calibrate the astrometry of each SW exposures. This information will be recorded in a table called tb_sw_astrometry and will be passed to the next step for LW grism astrometry correction.
Read LW direct imaging catalog, compare with SW extracted catalog for grism WCS correction¶
In this example, the astrometric catalog for IceAge data is the LW direct imaging catalog in F430M band. You can get access to it at http://gxn.as.arizona.edu/~sunfengwu/IceAge/mosaics/daofind_IceAge_NIRCam_F430M_source.cat.
In [35]:
'''Read astrometric catalog'''
tb_LW = ascii.read('/home/u24/fengwusun/IceAge/daofind_IceAge_NIRCam_F430M_source.cat')
tb_LW = tb_LW[tb_LW['flags'] == 0]
tb_LW[:3]
Out[35]:
In [37]:
'''*_cal.fits files in this band (F200W)'''
list_img_cal_this_band = list_rate_sw[list_sw_filter == 'F200W']
list_img_cal_this_band = np.array([calibrated_dir + '../direct/' + os.path.basename(x).replace('_rate.fits', '_cal.fits')
for x in list_img_cal_this_band])
In the following code block, you will run a loop to compute the x/y offset and rotation of four SW exposures associated with the LW grism exposure (e.g., A1-4 for module A) and save the astrometric shift information in the tb_sw_astrometry table. Note that module A/B might have different offsets so you have to group them separately.
In [38]:
'''SW astrometry table'''
tb_sw_astrometry = Table(names = ['expName','N_match', 'dRA', 'dDEC', 'dRA_err', 'dDEC_err', 'theta', 'theta_err'],
dtype = ['S40', 'i4', 'f4', 'f4' , 'f4', 'f4', 'f4', 'f4'])
tb_sw_astrometry.meta['comments'] = ['RA/Dec offsets are in arcsec and absolute',
'i.e., cos(Dec) projection has been considered.',
' >> dRA = RA_daofind - RA_gaia',
' >> dDEC = DEC_daofind - DEC_gaia',
'generated by F. Sun (steward) on %s' % time.strftime("%Y/%m/%d", time.localtime())]
for colname in tb_sw_astrometry.colnames[2:]: tb_sw_astrometry[colname].info.format = '.4f'
sigma_clip = SigmaClip(sigma = 2.)
'''Run the following loop to compute astrometric offsets for each group of SW exposures.'''
print('>>> Astrometry (RA/DEC) offsets calculation for %3d SW Frames' % len(list_img_cal_this_band))
for k in range(len(list_img_cal_this_band) // 4):
## For each cal.fits, get all the paths of csal.fits from all four detector.
tmp_img_cal_path = list_img_cal_this_band[k * 4]
tmp_img_cal_path_base = os.path.basename(tmp_img_cal_path)[:30]
tmp_obsnum = int(tmp_img_cal_path_base[7:10])
tmp_grism_LW_path = calibrated_dir + tmp_img_cal_path_base + 'long_rate_lv1.5.fits'
tmp_grism_hd_sci = fits.getheader(tmp_grism_LW_path, 1)
# if tmp_obsnum not in [125, 126]: continue
print('%3d >> %s ' % (k, tmp_img_cal_path_base))
tmp_img_cal_hd = fits.getheader(tmp_img_cal_path, 0)
tmp_detector = tmp_img_cal_hd['DETECTOR'].lower()
tmp_img_cal_sci_hd = fits.getheader(tmp_img_cal_path, 'sci')
tmp_time_obs = Time(tmp_img_cal_sci_hd['MJD-AVG'], format ='mjd')
tmp_detectors = np.array(['%s%d' % (tmp_detector[:-1], tmp_det) for tmp_det in [1, 2, 3, 4]])
tmp_img_cal_paths = np.array([tmp_img_cal_path.replace(tmp_detector, x) for x in tmp_detectors])
## Run DAOFIND to get a star catalog:
daofind = DAOStarFinder(fwhm = 3.0, threshold = 8.)
for j, tmp_img_cal_path in enumerate(tmp_img_cal_paths):
tmp_img_cal_sci_hd = fits.getheader(tmp_img_cal_path, 'sci')
'''Directly use prepared SW DAOFIND catalog'''
tmp_tb_daofind = ascii.read(calibrated_dir + 'astrom/%s' % os.path.basename(tmp_img_cal_path).replace('_cal.fits', '_daofind.cat').replace('_rate.fits', '_daofind.cat'),
format = 'ecsv')
'''no need to do daofind again'''
# tmp_wcs = wcs.WCS(tmp_img_cal_sci_hd)
# tmp_coord_center = SkyCoord(tmp_img_cal_sci_hd['crval1'], tmp_img_cal_sci_hd['crval2'], unit = (u.deg, u.deg))
# tmp_detect_img = np.nan_to_num(fits.getdata(tmp_img_cal_path, 'sci'))
# bkg_estimator = MedianBackground()
# bkg = Background2D(tmp_detect_img, (64, 64), filter_size = (5, 5),
# sigma_clip = SigmaClip(sigma = 3.), bkg_estimator = bkg_estimator)
# tmp_detect_std = sigma_clipped_stats(tmp_detect_img[100:500,100:500], sigma = 3, maxiters = 10)[-1]
# tmp_detect_img = (tmp_detect_img - bkg.background) / tmp_detect_std
# # tmp_detect_img[fits.getdata(tmp_img_cal_path, 'dq') >= 1] = 0
# # tmp_detect_img = np.nan_to_num(fits.getdata(tmp_img_cal_path, 'sci')/np.clip(fits.getdata(tmp_img_cal_path, 'err'), 1e-8, 1e8))
# tmp_tb_daofind = daofind(tmp_detect_img)
# tmp_tb_daofind['detector'] = tmp_detectors[j]
# tmp_tb_daofind['skycoord'] = wcs.utils.pixel_to_skycoord(tmp_tb_daofind['xcentroid'], tmp_tb_daofind['ycentroid'], tmp_wcs)
if j == 0 : tb_daofind = tmp_tb_daofind
else: tb_daofind = vstack((tb_daofind, tmp_tb_daofind))
tb_daofind = tb_daofind[np.argsort(tb_daofind['peak'])][-200:]
tmp_coord_daofind = tb_daofind['skycoord']
'''Astrometric Source RA & DECs'''
# tmp_coord_gaia = SkyCoord(tmp_RA, tmp_DEC, unit = (u.deg, u.deg))
# arg_LW_bright = np.where((tb_LW['MAG_AUTO'].data < 22.5) & (tb_LW['MAG_AUTO'].data > 16.5))[0]
# tmp_coord_LW = SkyCoord(tb_LW['X_WORLD'][arg_LW_bright], tb_LW['Y_WORLD'][arg_LW_bright],
# unit = (u.deg, u.deg))
arg_LW_bright = np.where((tb_LW['abmag'].data < 22.5) & (tb_LW['abmag'].data > 16.5))[0]
tmp_coord_LW = SkyCoord(tb_LW['RA'][arg_LW_bright], tb_LW['DEC'][arg_LW_bright],
unit = (u.deg, u.deg))
### rough center of the frame
tmp_coord_center = SkyCoord(np.median(tmp_coord_daofind.ra), np.median(tmp_coord_daofind.dec))
### only select sources close to the center
tmp_coord_LW = tmp_coord_LW[tmp_coord_center.separation(tmp_coord_LW) < 4 * u.arcmin]
## Cross match DAOFind Catalog with Astrometric Catalog
idx_daofind, d2d, _ = tmp_coord_LW.match_to_catalog_sky(tmp_coord_daofind)
idx_LW = np.where(d2d < 2.5 * u.arcsec)[0] ### Change the criterion if necessary
idx_daofind = idx_daofind[idx_LW]
## Compute RA/DEC offset
tmp_ra_offset = ((tmp_coord_daofind[idx_daofind].ra - tmp_coord_LW[idx_LW].ra)
* np.cos(tmp_coord_LW[idx_LW].dec)).to(u.arcsec).value
tmp_dec_offset = (tmp_coord_daofind[idx_daofind].dec - tmp_coord_LW[idx_LW].dec).to(u.arcsec).value
## sigma-clipped RA/DEC offset
tmp_ra_offset_clipped = sigma_clip(tmp_ra_offset)
tmp_dec_offset_clipped = sigma_clip(tmp_dec_offset)
arg_clipped = np.where((tmp_ra_offset_clipped.mask == False) & (tmp_dec_offset_clipped.mask == False))[0]
tmp_ra_offset_clipped, tmp_dec_offset_clipped = tmp_ra_offset[arg_clipped], tmp_dec_offset[arg_clipped]
tmp_ra_offset_med, tmp_ra_offset_std = np.median(tmp_ra_offset_clipped), np.std(tmp_ra_offset_clipped)
tmp_dec_offset_med, tmp_dec_offset_std = np.median(tmp_dec_offset_clipped), np.std(tmp_dec_offset_clipped)
## Compute rotation:
poly_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec, ydata = tmp_ra_offset)[0]
poly_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra, ydata = tmp_dec_offset)[0]
## iter-1
arg_dec_dRA = np.where(sigma_clip(np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value) - tmp_ra_offset).mask == False)[0]
poly_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec[arg_dec_dRA], ydata = tmp_ra_offset[arg_dec_dRA])[0]
arg_ra_dDEC = np.where(sigma_clip(np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value) - tmp_dec_offset).mask == False)[0]
poly_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra[arg_ra_dDEC], ydata = tmp_dec_offset[arg_ra_dDEC])[0]
## iter-2
arg_dec_dRA = np.where(sigma_clip(np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value) - tmp_ra_offset).mask == False)[0]
poly_dec_dRA, pcov_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec[arg_dec_dRA], ydata = tmp_ra_offset[arg_dec_dRA])
arg_ra_dDEC = np.where(sigma_clip(np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value) - tmp_dec_offset).mask == False)[0]
poly_ra_dDEC, pccov_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra[arg_ra_dDEC], ydata = tmp_dec_offset[arg_ra_dDEC])
perr_dec_dRA = np.sqrt(np.diag(pcov_dec_dRA))
perr_ra_dDEC = np.sqrt(np.diag(pccov_ra_dDEC))
## combine two direction
cos_factor = np.cos(np.deg2rad(np.median(tmp_coord_daofind[idx_daofind].dec.value)))
theta_as_deg = (poly_dec_dRA[0] / perr_dec_dRA[0] - poly_ra_dDEC[0] / perr_ra_dDEC[0]) / (1/perr_dec_dRA[0] + cos_factor / perr_ra_dDEC[0]) / 3600.
theta_as_deg_err = np.sum(perr_dec_dRA[0]**2 + (perr_ra_dDEC[0] / cos_factor)**2)**0.5 / 2. / 3600.
theta_as_deg = np.rad2deg(theta_as_deg)
theta_as_deg_err = np.rad2deg(theta_as_deg_err)
tmp_ra_offset_med = tmp_grism_hd_sci['CRVAL2'] * poly_dec_dRA[0] + poly_dec_dRA[1]
tmp_dec_offset_med = tmp_grism_hd_sci['CRVAL1'] * poly_ra_dDEC[0] + poly_ra_dDEC[1]
tmp_ra_offset_std = sigma_clipped_stats(tmp_ra_offset)[2]
tmp_dec_offset_std = sigma_clipped_stats(tmp_dec_offset)[2]
print(' RA_image - RA_ref = %.3f"±%.3f"' % (tmp_ra_offset_med, tmp_ra_offset_std / np.sqrt(len(tmp_ra_offset)) ))
print('DEC_image - DEC_ref = %.3f"±%.3f"' % (tmp_dec_offset_med, tmp_dec_offset_std / np.sqrt(len(tmp_ra_offset)) ))
print(' rotation theta = %.3f±%.3f deg' % (theta_as_deg, theta_as_deg_err))
tb_sw_astrometry.add_row([tmp_img_cal_path_base, len(arg_clipped),
tmp_ra_offset_med, tmp_dec_offset_med,
tmp_ra_offset_std, tmp_dec_offset_std,
theta_as_deg, theta_as_deg_err])
plt.close()
plt.subplots(1, 1, figsize = (5, 5))
plt.plot(tmp_ra_offset, tmp_dec_offset, marker = '.', color = 'gray', ls = 'none')
plt.plot(tmp_ra_offset_clipped, tmp_dec_offset_clipped, marker = '.', color = 'k', ls = 'none')
plt.plot(tmp_ra_offset_med, tmp_dec_offset_med, marker = '+', color = 'red', ms = 20, mew = 2, ls = 'none')
plt.axhline(0, color = 'grey', ls = '--')
plt.axvline(0, color = 'grey', ls = '--')
plt.gca().set(xlabel = r'$\Delta$RA', ylabel = r'$\Delta$Dec')
if np.max(np.abs(np.array([tmp_ra_offset_med, tmp_dec_offset_med]))) > 1:
plt.gca().set(aspect = 1, xlim = (-1.5, 1.5), ylim = (-1.5, 1.5),
xticks = np.arange(-1.5, 1.51, 0.5), yticks = np.arange(-1.5, 1.51, 0.3))
elif np.max(np.abs(np.array([tmp_ra_offset_med, tmp_dec_offset_med]))) > 0.5:
plt.gca().set(aspect = 1, xlim = (-1., 1.), ylim = (-1., 1.),
xticks = np.arange(-0.8, 1.1, 0.4), yticks = np.arange(-1., 1.1, 0.2))
else:
plt.gca().set(aspect = 1, xlim = (-0.5, 0.5), ylim = (-0.5, 0.5),
xticks = np.arange(-0.4, 0.51, 0.2), yticks = np.arange(-0.5, 0.51, 0.1))
corner_text(plt.gca(), s = ' RA_image - RA_ref = %.3f"±%.3f"' % (tmp_ra_offset_med, tmp_ra_offset_std / np.sqrt(len(tmp_ra_offset)) ),
loc = 1, fontsize = 12)
corner_text(plt.gca(), s = 'DEC_image - DEC_ref = %.3f"±%.3f"' % (tmp_dec_offset_med, tmp_dec_offset_std/ np.sqrt(len(tmp_ra_offset)) ),
loc = 4, fontsize = 12)
plt.title(tmp_img_cal_path_base, fontsize = 14, )
plt.tight_layout()
# plt.savefig('/home/u24/fengwusun/IceAge/SW_astrometry/%s_astrom.jpg' % tmp_img_cal_path_base, dpi = 100)
# break
In [ ]:
In [40]:
'''Checking the derived rotation, if necessary'''
fig, ax = plt.subplots(1, 2, figsize = (12, 6) )
ax[0].plot(tmp_coord_daofind[idx_daofind].dec, tmp_ra_offset, marker = '.', color = 'k', ls = 'none')
ax[0].set_ylim(-0.4, 0.4)
ax[1].plot(tmp_coord_daofind[idx_daofind].ra, tmp_dec_offset, marker = '.', color = 'k', ls = 'none')
ax[1].set_ylim(-0.4, 0.4)
poly_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec, ydata = tmp_ra_offset)[0]
poly_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra, ydata = tmp_dec_offset)[0]
## iter-1
arg_dec_dRA = np.where(sigma_clip(np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value) - tmp_ra_offset).mask == False)[0]
poly_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec[arg_dec_dRA], ydata = tmp_ra_offset[arg_dec_dRA])[0]
arg_ra_dDEC = np.where(sigma_clip(np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value) - tmp_dec_offset).mask == False)[0]
poly_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra[arg_ra_dDEC], ydata = tmp_dec_offset[arg_ra_dDEC])[0]
## iter-2
arg_dec_dRA = np.where(sigma_clip(np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value) - tmp_ra_offset).mask == False)[0]
poly_dec_dRA, pcov_dec_dRA = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].dec[arg_dec_dRA], ydata = tmp_ra_offset[arg_dec_dRA])
arg_ra_dDEC = np.where(sigma_clip(np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value) - tmp_dec_offset).mask == False)[0]
poly_ra_dDEC, pccov_ra_dDEC = optimize.curve_fit(f = linear, xdata = tmp_coord_daofind[idx_daofind].ra[arg_ra_dDEC], ydata = tmp_dec_offset[arg_ra_dDEC])
perr_dec_dRA = np.sqrt(np.diag(pcov_dec_dRA))
perr_ra_dDEC = np.sqrt(np.diag(pccov_ra_dDEC))
ax[0].plot(tmp_coord_daofind[idx_daofind].dec.value,
np.polyval(poly_dec_dRA, tmp_coord_daofind[idx_daofind].dec.value))
ax[1].plot(tmp_coord_daofind[idx_daofind].ra.value,
np.polyval(poly_ra_dDEC, tmp_coord_daofind[idx_daofind].ra.value))
ax[0].set(xlabel = 'Dec', ylabel = 'dRA')
ax[1].set(xlabel = 'RA', ylabel = 'dDEC')
plt.tight_layout()
print(poly_dec_dRA, poly_ra_dDEC)
Write astrometric catalog if necessary. The one that I got for Ice Age data (F200W image & F444W grism) is available at http://gxn.as.arizona.edu/~sunfengwu/IceAge/PID_1309_SW_wavecal_F200W_F430M.dat.
In [43]:
# ascii.write(tb_sw_astrometry, '/home/u24/fengwusun/IceAge/SW_astrometry/PID_1309_SW_wavecal_F200W_F430M.dat',
# overwrite = True)
In [42]:
# tb_sw_astrometry = ascii.read('/home/u24/fengwusun/IceAge/SW_astrometry/PID_1309_SW_wavecal_relative_astrom_obs24.dat')
print('0-16-50-84-100 percentiles of astrometric accuracy for RA direction [in mas]:')
print(np.round(np.percentile(tb_sw_astrometry['dRA_err'].data / np.sqrt(tb_sw_astrometry['N_match'].data - 1),
[0, 16, 50, 84, 100]) * 1e3, 2))
In [46]:
'''Evolution of astrometric offsets with exposures:'''
fig = plt.figure(constrained_layout = False, figsize = (12, 6))
ax_dict = fig.subplot_mosaic('''AACC
BBCC''')
ax_dict['A'].plot(tb_sw_astrometry['dRA']) # color = plt.cm.Reds(0.5))
ax_dict['B'].plot(tb_sw_astrometry['dDEC']) # color = plt.cm.Blues(0.5))
ax_dict['A'].set(ylabel = 'dRA [arcsec]', ylim = (-0.6, 0.6))
ax_dict['B'].set(ylabel = 'dDec [arcsec]', ylim = (-0.6, 0.6), xlabel = '$i_\mathrm{frame}$')
for i, tmp_dRA in enumerate(tb_sw_astrometry['dRA']):
if i%2 == 0: # mod-A
kwargs_dot = dict(color = plt.cm.Reds(i/len(tb_sw_astrometry)), marker = 'o', zorder = i)
else: # mod-B
kwargs_dot = dict(color = plt.cm.Blues(i/len(tb_sw_astrometry)), marker = 's', zorder = i)
ax_dict['C'].plot(tb_sw_astrometry['dRA'][i], tb_sw_astrometry['dDEC'][i], **kwargs_dot,
mec = 'silver', mew = 1, ms = 8)
ax_dict['C'].plot(tb_sw_astrometry['dRA'][::2], tb_sw_astrometry['dDEC'][::2],
ls = '-', color = plt.cm.Reds(0.5), zorder = 1)
ax_dict['C'].plot(tb_sw_astrometry['dRA'][1::2], tb_sw_astrometry['dDEC'][1::2],
ls = '-', color = plt.cm.Blues(0.5), zorder = 0)
ax_dict['C'].set(aspect = 1, ylabel = 'dDEC [arcsec]', xlabel = 'dRA [arcsec]')
corner_text(ax_dict['C'], s = 'modA', color = plt.cm.Reds(0.5), loc = 2, weight = 'semibold')
corner_text(ax_dict['C'], s = 'modB', color = plt.cm.Blues(0.5), loc = 1, weight = 'semibold')
fig.text(0.5, 0.98, 'PID-1309 SW Relative Astrometric Offset', ha = 'center', va = 'top', fontsize = 18)
plt.tight_layout()
plt.subplots_adjust(top = 0.93)
# fig.savefig('/home/u24/fengwusun/IceAge/PID_1309_SW_relative_astrom_F200W.pdf', dpi = 100)
In [47]:
'''Astrometric offsets between module A and B'''
plt.errorbar(np.arange(len(tb_sw_astrometry)/2),
((tb_sw_astrometry['dRA'][::2] - tb_sw_astrometry['dRA'][1::2])**2 +
(tb_sw_astrometry['dDEC'][::2] - tb_sw_astrometry['dDEC'][1::2])**2)**0.5,
yerr = tb_sw_astrometry['dDEC_err'][::2] * np.sqrt(2) / np.sqrt(tb_sw_astrometry['N_match'][::2] - 1),
ls = '-', zorder = 0)
plt.gca().set(xlabel = 'Exp NUM', ylabel = 'modA - modB OFFSET', ylim = (-0.2, 0.5))
Out[47]:
Check source position predicted by fits WCS + Measured Astrometry on SW images¶
This method will also be used to predict source positions on grism data
In [50]:
## Select one SW exposure to look at:
k = 0
tmp_img_cal_path = list_img_cal_this_band[k * 4 + 1]
tmp_img_cal_path_base = os.path.basename(tmp_img_cal_path)[:30]
## astrometric info:
item_sw_astrom = tb_sw_astrometry[tb_sw_astrometry['expName'] == tmp_img_cal_path_base][0]
print('%3d >> %s ' % (k, tmp_img_cal_path))
tmp_img_cal_hd = fits.getheader(tmp_img_cal_path, 0)
tmp_detector = tmp_img_cal_hd['DETECTOR'].lower()
tmp_img_cal_sci_hd = fits.getheader(tmp_img_cal_path, 'sci')
tmp_img_cal_wcs = wcs.WCS(tmp_img_cal_sci_hd)
time_obs = Time(tmp_img_cal_sci_hd['MJD-AVG'], format ='mjd')
tmp_detectors = np.array(['%s%d' % (tmp_detector[:-1], tmp_det) for tmp_det in [1, 2, 3, 4]])
tmp_img_cal_paths = np.array([tmp_img_cal_path.replace(tmp_detector, x) for x in tmp_detectors])
'''astrometric catalog loaded'''
# tmp_RA, tmp_DEC = tb_LW['ALPHA_J2000'].data, tb_LW['DELTA_J2000'].data
tmp_RA, tmp_DEC = tb_LW['RA'].data, tb_LW['DEC'].data
# tmp_RA_astrom = tmp_RA + item_sw_astrom['dRA'] / np.cos(np.deg2rad(tmp_DEC)) / 3600.
# tmp_DEC_astrom = tmp_DEC + item_sw_astrom['dDEC'] / 3600.
# tmp_img_cal_sci_hd['crval1'] -= item_sw_astrom['dRA'] / np.cos(np.deg2rad(tmp_img_cal_sci_hd['crval2'])) / 3600.
# tmp_img_cal_sci_hd['crval2'] -= item_sw_astrom['dDEC'] / 3600.
'''astrometric info applied on WCS:'''
tmp_cd_matrix = np.array([[tmp_img_cal_sci_hd['CD1_1'], tmp_img_cal_sci_hd['CD1_2']],
[tmp_img_cal_sci_hd['CD2_1'], tmp_img_cal_sci_hd['CD2_2']]])
tmp_rot_matrix = np.array([[np.cos(np.deg2rad(-item_sw_astrom['theta'])), - np.sin(np.deg2rad(-item_sw_astrom['theta']))],
[np.sin(np.deg2rad(-item_sw_astrom['theta'])), np.cos(np.deg2rad(-item_sw_astrom['theta']))]])
tmp_cd_matrix = np.matmul(tmp_cd_matrix, tmp_rot_matrix)
tmp_img_cal_sci_hd['CD1_1'] = tmp_cd_matrix[0][0]
tmp_img_cal_sci_hd['CD1_2'] = tmp_cd_matrix[0][1]
tmp_img_cal_sci_hd['CD2_1'] = tmp_cd_matrix[1][0]
tmp_img_cal_sci_hd['CD2_2'] = tmp_cd_matrix[1][1]
# ## Read latest SiAF for NIRCam
# nircam_siaf = pysiaf.Siaf('NIRCam')
# tmp_siaf = nircam_siaf['%s_%s' % (tmp_img_cal_hd['detector'], tmp_img_cal_hd['SUBARRAY'])]
# try: rotation = item_sw_astrom['theta']
# except KeyError: rotation = 0
'''
Convert coordinate of astrometric sources to pixel_x/y based on WCS
'''
plt.figure(figsize = (12, 12))
tmp_sci_img = fits.getdata(tmp_img_cal_path, 'sci')
plt.imshow(tmp_sci_img, vmin = np.nanpercentile(tmp_sci_img, 40), vmax = np.nanpercentile(tmp_sci_img, 98), origin = 'lower')
tmp_wcs = wcs.WCS(tmp_img_cal_sci_hd)
pixelx, pixely = wcs.utils.skycoord_to_pixel(SkyCoord(tmp_RA, tmp_DEC, unit = (u.deg, u.deg)), tmp_wcs)
arg_gaia_plot = np.where((pixelx > 0) & (pixelx < 2048) & (pixely > 0) & (pixely < 2048))[0]
plt.plot(pixelx[arg_gaia_plot], pixely[arg_gaia_plot], marker = 'o', mfc = 'none', mec = 'salmon', ms = 10, ls = 'none', mew = 2)
# tmp_asdf_distortion = asdf.open(crds.getreferences(tmp_img_cal_hd,
# reftypes = ['distortion'], ignore_cache = False)['distortion'])
# loc_v2, loc_v3 = tmp_asdf_distortion['model'](pixelx, pixely)
# pixelx, pixely = tmp_siaf.tel_to_sci(loc_v2, loc_v3)
# plt.plot(pixelx[arg_gaia_plot], pixely[arg_gaia_plot], marker = 'o', mfc = 'none', mec = 'w', ms = 10, ls = 'none', mew = 2)
# plt.xlim(200, 700);plt.ylim(200, 700)
# plt.xlim(1300, 1800);plt.ylim(1300, 1800)
Out[50]:
In [ ]:
5. Genertae Source Catalog for each LW grism exposure¶
With astrometric information in hand, we can start to generate the catalog that include sources of our interest for each LW grism exposures. With these catalogs (I name them as "POM-applied catalogs" because POM selection effects are applied), you can then run the code blocks to extract 2D spectra from individual LW grism frames, co-add them together to produce the final 2D grism spectra.
First of all, let us grab the reduced LW grism images (_rate_lv1.5.fits) again:
In [51]:
'''LW Grism images and corresponding astrometric information:'''
# all_rate_list = glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/grism/jw0130902[3-4]*_rate_lv1.5.fits')
# tb_sw_astrometry = ascii.read('/home/u24/fengwusun/IceAge/SW_astrometry/PID_1309_SW_wavecal_F150W_F410M.dat')
all_rate_list = glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/grism/jw0130902[5-6]*_rate_lv1.5.fits')
tb_sw_astrometry = ascii.read('/home/u24/fengwusun/IceAge/SW_astrometry/PID_1309_SW_wavecal_F200W_F430M.dat')
all_rate_list.sort()
all_rate_list = np.array(all_rate_list)
print(len(all_rate_list), 'wcs corrected fits files')
'''read filter, module and pupil'''
all_filters, all_module, all_pupil = [], [], []
all_roll_angle, all_target = [], []
for i, x in enumerate(all_rate_list):
tmp_header = fits.getheader(x)
all_filters.append(tmp_header['filter'])
all_module.append(tmp_header['module'])
all_pupil.append(tmp_header['pupil'])
all_roll_angle.append(fits.getheader(x, 1)['ROLL_REF'])
all_target.append(tmp_header['TARGPROP'])
all_filters, all_module, all_pupil = np.array(all_filters), np.array(all_module), np.array(all_pupil)
all_roll_angle, all_target = np.array(all_roll_angle), np.array(all_target)
print(np.unique(all_filters), np.unique(all_module), np.unique(all_target))
### change target prefix if you like:
tmp_filter = np.unique(all_filters)[0]
if tmp_filter == 'F322W2': tmp_target = 'CHAMMS1-H2O'
if tmp_filter == 'F444W': tmp_target = 'CHAMMS1-CO'
Ice Age data consist of F322W2 and F444W grism spectra taken with both dispersion direction (R and C). In the following example, we input the source catalog extracted in both F410M (imaged after F322W2 grism exposure) and F430M (imaged after F444W grism exposure) band as master catalog for spectral extraction.
In [53]:
'''Input master catalog for spectral extraction:'''
tb_f322w2 = ascii.read('/home/u24/fengwusun/IceAge/daofind_IceAge_NIRCam_F410M_source.cat')
tb_f444w = ascii.read('/home/u24/fengwusun/IceAge/daofind_IceAge_NIRCam_F430M_source.cat')
tb_f322w2 = tb_f322w2[tb_f322w2['flags'] == 0]
tb_f444w = tb_f444w[tb_f444w['flags'] == 0]
'''If necessary, select source based on morphology and brigtness'''
# tb_f322w2 = tb_f322w2[(tb_f322w2['MAG_AUTO'] < 18.5) & (tb_f322w2['CLASS_STAR'] > 0.5)]
### read coordniates
coord_f322w2 = SkyCoord(tb_f322w2['RA'], tb_f322w2['DEC'], unit = (u.deg, u.deg))
coord_f444w = SkyCoord(tb_f444w['RA'], tb_f444w['DEC'], unit = (u.deg, u.deg))
print(len(coord_f322w2), len(coord_f444w))
In [55]:
'''Columnames of my catalog'''
print(tb_f444w.colnames)
tb_f444w[:4]['id'].data
Out[55]:
For each exposure, apply astrometric correction and generate source catalog for spectral extraction¶
In this step, there is one tricky thing called filter offset. As you might have recalled, these offsets are introduced by different SW+LW filter combination. In commissioning, we used F212N + F444W for F444W grism spectroscopy, but Ice Age used F200W + F444W. Therefore, the spectral tracing model derived in commissioning could be offset from the actual ones in IceAge data due to the difference in SW filters. This should be a global shift (dx, dy) which is different for different modules but consistent for all exposures. I believe that this has not been well calibrated so far, so you will need to enter some dx, dy value manually. You can use the script a few blocks below the following code block to check the filter offsets in current CRDS, but my general advice is to change these dx, dy number manually after you successfully extract some 2D spectra and measure the real offsets.
In [60]:
'''POM applied catalogs'''
### `all_rate_list`: a list of paths to the grism exposure fits files
for i, tmp_rate_path in enumerate(all_rate_list[:]):
'''
## The following commands are used to read the astrometric offsets which is measured with SW detectors.
## In real life, the pointing accuracy of the telescope may not be stable and the absolute RA and DEC
## offset can be floating. We will need to measure the astrometric shift from the SW exposures associated
## with the grism exposures, and correct this factor for the RA and DEC in each source catalog associated
## with the grism exposures.
'''
tmp_rate_path_base = tmp_rate_path.split('/')[-1].split('long_rate')[0]
print('[%3d]' % i, tmp_rate_path_base)
tmp_grism_hd_1st = fits.getheader(tmp_rate_path, 0)
tmp_grism_hd_sci = fits.getheader(tmp_rate_path, 'sci')
tmp_filter = tmp_grism_hd_1st['filter']
tmp_loc_v2, tmp_loc_v3 = tmp_grism_hd_sci['V2_REF'], tmp_grism_hd_sci['V3_REF']
### for different grism filters, read different catalog
if tmp_filter == 'F444W': list_coords, tb_sex = coord_f444w, tb_f444w
elif tmp_filter == 'F322W2': list_coords, tb_sex = coord_f322w2, tb_f322w2
# if tmp_grism_hd_1st['filter'] != 'F444W': continue
'''Accurate RA & Dec of Direct-Imaging sources '''
tmp_RA, tmp_DEC = list_coords.ra.to(u.deg).value, list_coords.dec.to(u.deg).value
### if the catalog is gaia star, considering proper motion:
# time_obs = Time(tmp_grism_hd_sci['MJD-AVG'], format ='mjd')
# dRA = ((time_obs - time_gaia).to(u.yr).value * tb_gaia['pmRA'].data * u.mas / np.cos(np.deg2rad(tb_gaia['DE_ICRS']))).to(u.deg).value
# dDec = ((time_obs - time_gaia).to(u.yr).value * tb_gaia['pmDE'].data * u.mas).to(u.deg).value
# tmp_RA, tmp_DEC = tb_gaia['RA_ICRS'].data.data + dRA, tb_gaia['DE_ICRS'].data.data + dDec
'''Apply astrometric error measured on SW images'''
item_sw_astrom = tb_sw_astrometry[tb_sw_astrometry['expName'] == tmp_rate_path_base]
if len(item_sw_astrom) > 0:
## if there is astrometric information, consider that
item_sw_astrom = item_sw_astrom[0]
tmp_grism_hd_sci['CRVAL1'] -= item_sw_astrom['dRA'] / np.cos(np.deg2rad(tmp_grism_hd_sci['CRVAL2'])) / 3600.
tmp_grism_hd_sci['CRVAL2'] -= item_sw_astrom['dDEC'] / 3600.
tmp_cd_matrix = np.array([[tmp_grism_hd_sci['CD1_1'], tmp_grism_hd_sci['CD1_2']],
[tmp_grism_hd_sci['CD2_1'], tmp_grism_hd_sci['CD2_2']]])
tmp_rot_matrix = np.array([[np.cos(np.deg2rad(-item_sw_astrom['theta'])), - np.sin(np.deg2rad(-item_sw_astrom['theta']))],
[np.sin(np.deg2rad(-item_sw_astrom['theta'])), np.cos(np.deg2rad(-item_sw_astrom['theta']))]])
tmp_cd_matrix = np.matmul(tmp_cd_matrix, tmp_rot_matrix)
tmp_grism_hd_sci['CD1_1'] = tmp_cd_matrix[0][0]
tmp_grism_hd_sci['CD1_2'] = tmp_cd_matrix[0][1]
tmp_grism_hd_sci['CD2_1'] = tmp_cd_matrix[1][0]
tmp_grism_hd_sci['CD2_2'] = tmp_cd_matrix[1][1]
tmp_RA_astrom = tmp_RA # + item_sw_astrom['dRA'] / np.cos(np.deg2rad(tmp_DEC)) / 3600.
tmp_DEC_astrom = tmp_DEC # + item_sw_astrom['dDEC'] / 3600.
else:
## if no astrometric information avaialble, skip this step and directly use RA and DEC
print('Warning: No SW images found for %s, skip astrom correction! ' % tmp_rate_path_base)
tmp_RA_astrom = tmp_RA
tmp_DEC_astrom = tmp_DEC
try: rotation = item_sw_astrom['theta']
except KeyError: rotation = 0
'''RA/DEC to grism-frame x and y'''
tmp_grism_wcs = wcs.WCS(tmp_grism_hd_sci)
idx_this_field = np.where((np.abs((tmp_RA_astrom - tmp_grism_hd_sci['crval1']) * np.cos(np.deg2rad(tmp_grism_hd_sci['crval2']))) < 3 / 60.) &
(np.abs(tmp_DEC_astrom - tmp_grism_hd_sci['crval2']) < 3 / 60.))[0]
pixelx, pixely = wcs.utils.skycoord_to_pixel(SkyCoord(tmp_RA_astrom[idx_this_field],
tmp_DEC_astrom[idx_this_field], unit = (u.deg, u.deg)),
tmp_grism_wcs)
'''This is the place to add small offsets from all kind of reasons (e.g., filter offset)'''
if tmp_grism_hd_1st['module'] == 'A' : dx, dy = 0, 0 ### +2.0, -2.5
else: dx, dy = 0, 0
pixelx, pixely = pixelx + dx, pixely + dy
tb_sex = tb_sex[idx_this_field]
'''for sextractor-based catalog '''
# tmp_mag_auto, tmp_magerr_auto = tb_sex['MAG_AUTO'].data, tb_sex['MAGERR_AUTO'].data
# tb_pom_applied = Table(data = [tb_sex['NUMBER'].data,
# tmp_RA[idx_this_field], tmp_DEC[idx_this_field],
# pixelx, pixely, tmp_mag_auto, tmp_magerr_auto,
# tb_sex['A_IMAGE'], tb_sex['B_IMAGE'], tb_sex['THETA_IMAGE'],
# tb_sex['FLAGS'], tb_sex['CLASS_STAR']],
# names = ['Index', 'ra', 'dec', 'pixel_x', 'pixel_y', 'MAG_AUTO', 'MAGERR_AUTO',
# 'A_IMAGE', 'B_IMAGE', 'THETA_IMAGE', 'FLAGS', 'CLASS_STAR'])
'''for DAOFIND-based catalog '''
tmp_mag_auto, tmp_magerr_auto = tb_sex['abmag'].data, tb_sex['abmag_err'].data
tb_pom_applied = Table(data = [tb_sex['id'].data,
tmp_RA[idx_this_field], tmp_DEC[idx_this_field],
pixelx, pixely, tmp_mag_auto, tmp_magerr_auto, tb_sex['flags']],
names = ['Index', 'ra', 'dec', 'pixel_x', 'pixel_y', 'MAG_AUTO', 'MAGERR_AUTO', 'FLAGS'])
# arg_is_pickoff = np.where(is_pickoff(pixelx, pixely, module = tmp_grism_hd_1st['module']))[0]
arg_is_pickoff = np.where(is_pickoff_PS(pixelx, pixely, module = tmp_grism_hd_1st['module'],
filter = tmp_grism_hd_1st['filter'], pupil = tmp_grism_hd_1st['pupil'][-1]) > 0.05)[0]
tb_pom_applied = tb_pom_applied[arg_is_pickoff]
tb_pom_applied['ra'].info.format = '.6f'
tb_pom_applied['dec'].info.format = '.6f'
tb_pom_applied['pixel_x'].info.format = '.3f'
tb_pom_applied['pixel_y'].info.format = '.3f'
tb_pom_applied['MAG_AUTO'].info.format = '.3f'
tb_pom_applied['MAGERR_AUTO'].info.format = '.3f'
try:
tb_pom_applied['A_IMAGE'].info.format = '.3f'
tb_pom_applied['B_IMAGE'].info.format = '.3f'
tb_pom_applied['THETA_IMAGE'].info.format = '.2f'
tb_pom_applied['CLASS_STAR'].info.format = '.3f'
except KeyError: pass
'''Change the path of catalogs to be saved:'''
path_tb_pom_applied = '/home/u24/fengwusun/IceAge/POM_catalog/' + tmp_rate_path_base + '_%s_dirimg_sources.list' % tmp_filter
ascii.write(tb_pom_applied, path_tb_pom_applied, overwrite = True)
Optional: Check filter offset:¶
This should have been corrected in stage-2 using STScI CRDS files. However, I stil see some offsets in modA, weird!
In [57]:
### select a reduced grism exposure:
tmp_rate_path = '/xdisk/egami/fengwusun/ERSs/IceAge/grism/jw01309025002_02101_00012_nrcblong_rate_lv1.5.fits'
tmp_rate_path_base = tmp_rate_path.split('/')[-1].split('long_rate')[0]
### SW header
tmp_sw_cal_path = os.path.dirname(tmp_rate_path) + '/../direct/' + tmp_rate_path_base + '1_cal.fits'
tmp_sw_cal_hd = fits.getheader(tmp_sw_cal_path)
### LW header
tmp_grism_hd_1st = fits.getheader(tmp_rate_path, 0)
### offsets in LW & SW
tmp_tb_foff_lw = Table(asdf.open(crds.getreferences(tmp_grism_hd_1st)['filteroffset'])['filters'])
tmp_tb_foff_sw = Table(asdf.open(crds.getreferences(tmp_sw_cal_hd)['filteroffset'])['filters'])
### Filter offset of F212N, which is used in Commissioning
item_F212N = tmp_tb_foff_sw[(tmp_tb_foff_sw['filter'] == 'F212N') & (tmp_tb_foff_sw['pupil'] == tmp_sw_cal_hd['pupil'])][0]
### Filter offset of the filter used in the current SW exposures
item_sw = tmp_tb_foff_sw[(tmp_tb_foff_sw['filter'] == tmp_sw_cal_hd['filter']) & (tmp_tb_foff_sw['pupil'] == tmp_sw_cal_hd['pupil'])][0]
In [59]:
item_F212N, item_sw
# tmp_tb_foff_sw[(tmp_tb_foff_sw['filter'] == 'F200W') & (tmp_tb_foff_sw['pupil'] == tmp_sw_cal_hd['pupil'])][0]
Out[59]:
In [739]:
'''Sometimes LW filters also introduce offset (e.g., F356W - F322W2). Check that if necessary'''
# tmp_tb_foff_lw[(tmp_tb_foff_lw['filter'] == 'F356W') & (tmp_tb_foff_lw['pupil'] == 'CLEAR')]
6. 2D Spectral Extraction¶
The following step is to extract grism spectra with the POM-applied catalogs and reduced grism data you produced above. At this step, you will need to load the following functions:
- Spectral Tracing Model
Now we assume $(x_0, y_0)$ are direct imaging source position, $(x_s, y_s)$ are spectral pixel position at wavelength $\lambda_s$, and $dx = x_s - x_0$, $dy = y_s - y_0$. The spectral tracing model gives $dy(x_0, y_0, dx)$, i.e., the spectral tracing. The following plot shows the spectral tracing that we derived for F322W2 module A grism R using the commissioning data taken for the LMC field. Red squares are direct imaging positions of sources taken in multiple exposures, and black points are the corresponding spectral traces (curvature is exaggerated). The plot of dx versus dy, color-coded by the source position, is shown on the right. Note the spectra with $dx \sim 1000$ are second-order, which have not been fully characterized at this step.
The spectral tracing function is defined as fit_disp_order23() (see the beginning of this notebook). The best-fit parameters for these functions in different filter/module/grism can be found at http://u.arizona.edu/~fengwusun/data/dy_dx_tracing_model/. For filters at 2.5-4.0 µm (e.g., F356W), please simply use the models for F322W2; for filters at 4.0-5.0µm (e.g., F410M), please use the models for F444W.
These models have a typical accuracy of RMS<0.1 pixel if astrometric and filter offsets are corrected. Below is a plot showing the $dy(dx=0)$ map for F444W module A grism R. The dy of my model is color-coded as the background image, and squares with filled colors is what I directly measured from the spectral traces in the LMC field. You can see a field variation of $dy(dx=0)$ of 1-2 pixel, which has been well modeled with my spectral tracing function.
- Grism Dispersion model
The grism dispersion models give the relation of $dx(x_0, y_0, \lambda_s)$, i.e., the dispersion at given wavelength for source at $(x_0, y_0)$. In the commissioning, this was characterized by the F322W2 and F444W spectra of post-AGB star IRAS 05248-7007 in the LMC field.
The grism dispersion function is defined as fit_disp_order32() (or func_fit_wave(); see the beginning of this notebook). The best-fit parameters for these functions in different filter/module/grism can be found at http://u.arizona.edu/~fengwusun/data/dx_wave_dispersion_model/. The function for the same grism + module combination applies for all filters.
The accuracy of wavelength calibration is on the order of RMS=0.1-0.2 pixel. However, I would like to note that the current wavelength calibration data is still quite limited (~200 data points fit for a 16-parameter functions), and I would expect the accuracy gets further improved with cycle-1 calibration program.
- Sensitivity Function
All of the grism transmission curves that I have produced are available at http://u.arizona.edu/~fengwusun/data/sensitivity_model/. These files were derived using the Cycle-1 WFSS flux cal observations obtained in September, 2022. All of the grism R flux calibration is based on a G dwarf (P330-E; PID: 1538), a white dwarf (G191-B2B; PID: 1537) and an A-type star (J1743045; PID: 1536) with six exposures in total. The C grism observation of the white dwarf is corrupted by nearby stars, so most of the flux calibration of grism C is only based on P330E and J1743045 (4 individual integrations in total). The derived grism sensitivity functions in all filter+grism+module combinations are shown as follows:
The three columns of each grism sensitivity table are wavelength (in µm), sensitivity (in DN/s/Jy/pixel) and uncertainty of sensitivity. The sensitivity is the total throughput (OTE QE filter * grism, etc.) and sums the counts across the perpendicular (i.e. spatial) direction of the spectral trace, so it is in unit of DN/s/Jy/pixel. The dispersion is around 0.97 nm/pixel, so it is possible to convert the sensitivity to those in unit like (DN/s)/(erg/s/cm^2/Å)/nm. Note that the dispersion (nm/pixel or Å/pixel) is a function of wavelength and direct imaging position $(x0, y0)$, but 0.97 nm/pixel is a good characteristic dispersion.
Load the grism tracing/dispersion/sensitivity functions and POM-applied catalogs:¶
In [63]:
### select the filter that we will work on:
tmp_filter = 'F444W'
if tmp_filter == 'F444W': list_coords, tb_sex = coord_f444w, tb_f444w
if tmp_filter == 'F322W2': list_coords, tb_sex = coord_f322w2, tb_f322w2
### Interested wavelength Range
if tmp_filter == 'F444W': WRANGE = np.array([3.8, 5.1])
elif tmp_filter == 'F322W2': WRANGE = np.array([2.4, 4.1])
elif tmp_filter == 'F356W': WRANGE = np.array([3.1, 4.0])
elif tmp_filter == 'F277W': WRANGE = np.array([2.4, 3.1])
### Spectral tracing parameters:
if tmp_filter in ['F277W', 'F356W']: disp_filter = 'F322W2'
else: disp_filter = tmp_filter
tb_order23_fit_AR = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, 'A', 'R'))
fit_opt_fit_AR, fit_err_fit_AR = tb_order23_fit_AR['col0'].data, tb_order23_fit_AR['col1'].data
tb_order23_fit_BR = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, 'B', 'R'))
fit_opt_fit_BR, fit_err_fit_BR = tb_order23_fit_BR['col0'].data, tb_order23_fit_BR['col1'].data
tb_order23_fit_AC = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, 'A', 'C'))
fit_opt_fit_AC, fit_err_fit_AC = tb_order23_fit_AC['col0'].data, tb_order23_fit_AC['col1'].data
tb_order23_fit_BC = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, 'B', 'C'))
fit_opt_fit_BC, fit_err_fit_BC = tb_order23_fit_BC['col0'].data, tb_order23_fit_BC['col1'].data
### grism dispersion parameters:
tb_fit_displ_AR = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % ('A', "R"))
w_opt_AR, w_err_AR = tb_fit_displ_AR['col0'].data, tb_fit_displ_AR['col1'].data
tb_fit_displ_BR = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % ('B', "R"))
w_opt_BR, w_err_BR = tb_fit_displ_BR['col0'].data, tb_fit_displ_BR['col1'].data
tb_fit_displ_AC = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % ('A', "C"))
w_opt_AC, w_err_AC = tb_fit_displ_AC['col0'].data, tb_fit_displ_AC['col1'].data
tb_fit_displ_BC = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % ('B', "C"))
w_opt_BC, w_err_BC = tb_fit_displ_BC['col0'].data, tb_fit_displ_BC['col1'].data
### list of module/pupil and corresponding tracing/dispersion function:
list_mod_pupil = np.array(['AR', 'BR', 'AC', 'BC'])
list_fit_opt_fit = np.array([fit_opt_fit_AR, fit_opt_fit_BR, fit_opt_fit_AC, fit_opt_fit_BC])
list_w_opt = np.array([w_opt_AR, w_opt_BR, w_opt_AC, w_opt_BC])
### Sensitivity curve:
dir_fluxcal = '/home/u24/fengwusun/Comissioning/cycle1_cal_program/fluxcal/product/'
tb_sens_AR = ascii.read(dir_fluxcal + '%s_mod%s_grism%s_sensitivity.dat' % (disp_filter, 'A', 'R'))
tb_sens_BR = ascii.read(dir_fluxcal + '%s_mod%s_grism%s_sensitivity.dat'% (disp_filter, 'B', 'R'))
tb_sens_AC = ascii.read(dir_fluxcal + '%s_mod%s_grism%s_sensitivity.dat' % (disp_filter, 'A', 'C'))
tb_sens_BC = ascii.read(dir_fluxcal + '%s_mod%s_grism%s_sensitivity.dat'% (disp_filter, 'B', 'C'))
f_sens_AR = interpolate.UnivariateSpline(tb_sens_AR['wavelength'], tb_sens_AR['DN/s/Jy'], ext = 'zeros', k = 1, s = 1e2)
f_sens_BR = interpolate.UnivariateSpline(tb_sens_BR['wavelength'], tb_sens_BR['DN/s/Jy'], ext = 'zeros', k = 1, s = 1e2)
f_sens_AC = interpolate.UnivariateSpline(tb_sens_AC['wavelength'], tb_sens_AC['DN/s/Jy'], ext = 'zeros', k = 1, s = 1e2)
f_sens_BC = interpolate.UnivariateSpline(tb_sens_BC['wavelength'], tb_sens_BC['DN/s/Jy'], ext = 'zeros', k = 1, s = 1e2)
list_f_sens = (f_sens_AR, f_sens_BR, f_sens_AC, f_sens_BC)
'''list of source list (POM-applied catalog)'''
list_sl = glob.glob('/home/u24/fengwusun/IceAge/POM_catalog/jw01309*_nrc[a-b]_%s_dirimg_sources.list' % tmp_filter)
list_sl.sort()
list_rate_this = np.array([os.path.dirname(all_rate_list[0]) + '/' + os.path.basename(list_sl[i]).split('_%s' % tmp_filter)[0] + 'long_rate_lv1.5.fits'
for i in range(len(list_sl))])
### apply magnitude cut to the input catalog?
mag_cut = 18
Index_interest = tb_sex[tb_sex['abmag'] <= mag_cut]['id'].data
# Index_interest = tb_sex[tb_sex['MAG_AUTO'] <= mag_cut]['NUMBER'].data
# Index_interest = tb_sex[tb_sex['FLUX_AUTO'] >= 10**(-0.4 * (mag_cut - nircam_LW_ZP))]['NUMBER'].data
print(len(Index_interest), 'sources selected')
list_sl_table = [ascii.read(x) for x in list_sl]
'''
Generate two dictionaries which record the grism frame ID and path for each source to be extracted
'''
grism_frame_ID_per_source = dict() # dictionary of the ID of exposure (0-95) for each source (0-5000+)
grism_frame_path_per_source = dict() # dictionary of the path of exposure (*/jw*.fits) for each source (0-5000+)
for idx in Index_interest[:]:
grism_frame_ID_per_source['%s' % idx] = []
grism_frame_path_per_source['path_%s' % idx] = []
for i, tmp_sl_table in enumerate(list_sl_table):
for idx in tmp_sl_table['Index'].data:
try:
grism_frame_ID_per_source['%s' % idx].append(i)
grism_frame_path_per_source['path_%s' % idx].append(list_rate_this[i])
except KeyError: continue
N_frame_per_source = np.array([len(grism_frame_ID_per_source[x]) for x in grism_frame_ID_per_source.keys()])
The following functions are used to extract 2d spectra, resample / coadd 2d spectra and save them.¶
In [66]:
def extract_2d_spec(img, WRANGE, x0, y0, dxs, dys, wave, img_wht, img_dq, header,
aper = 10., pupil = 'R', n_procs = 4):
'''
Extract 2D spectrum from the Grism Image
-----------------------------------------------
Parameters
----------
img : `~numpy.ndarray`
2D dispersed slitless spectroscopic image
GConf: `~grismconf.grismconf.Config`
Configuration file of grism (taken from `grismconf`)
x0, y0 : float
Reference position (i.e., in direct image)
dxs, dys : `~numpy.ndarray`
First derivative of the trace x/y-coordinates with respect to 't' (taken from `grismconf`)
wave:
Wavelength of extracted spectrum
img_wht : `~numpy.ndarray`
Weighting image (1/rms^2) of `img`
img_dq : `~numpy.ndarray`
DQ (data-quality) image (1/rms^2) of `img`
header: `fits.Header()` object
Header of 2D dispersed slitless spectroscopic image
aper: float
Aperture radius for 2D spectra extraction (unit: arcsec, default: 1.25)
pix_scale: float
Pixel scale in unit of arcsec. If None then computed from header
pupil: 'R' or 'C'
pupil of grism ('R' or 'C')
Returns
-------
hdul : `~astropy.io.fits.HDUList()` object
The fits files of extracted 2D spectra, including a primary frame, a sci frame,
(scientific), a wht frame (weighting) and a dq frame (data quality).
'''
w_min, w_max = WRANGE
x_on_G_img = dxs + x0
y_on_G_img = dys + y0
if pupil not in ['C', 'R']: raise KeyError('pupil of grism should be either "R" or "C"! ')
aper_int_pix = int(aper)
### region that we can extract spectra:
if pupil == 'R':
args_eff = np.where((wave >= w_min) & (wave <= w_max)
& (x_on_G_img >= 5) & (x_on_G_img <= 2047 - 6)
& (y_on_G_img >= aper_int_pix) & (y_on_G_img <= 2047 - 6 - aper_int_pix) )
elif pupil == 'C':
args_eff = np.where((wave >= w_min) & (wave <= w_max)
& (x_on_G_img >= aper_int_pix) & (x_on_G_img <= 2047 - 6 - aper_int_pix)
& (y_on_G_img >= 5) & (y_on_G_img <= 2047 - 6) )
if np.size(args_eff) <= 20: ## spectrum is too short
raise ValueError('No spectrum can be extracted from the 2D Image')
tmp_spec_2d = np.zeros((len(args_eff[0]), aper_int_pix*2+1))
tmp_wht_2d = np.zeros((len(args_eff[0]), aper_int_pix*2+1))
tmp_dq_2d = np.zeros((len(args_eff[0]), aper_int_pix*2+1))
for i, j in enumerate(args_eff[0]):
if pupil == 'R':
tmp_x, tmp_y1, tmp_y2 = int(x_on_G_img[j]), int(y_on_G_img[j] - aper_int_pix - 1), int(y_on_G_img[j] + aper_int_pix + 2)
img.T[tmp_x, tmp_y1:tmp_y2] = ndimage.shift(img.T[tmp_x, tmp_y1:tmp_y2], -(y_on_G_img[j]%1), order = 1, mode ='wrap')
img_wht.T[tmp_x, tmp_y1:tmp_y2] = ndimage.shift(img_wht.T[tmp_x, tmp_y1:tmp_y2], -(y_on_G_img[j]%1), order = 1, mode ='wrap')
img_dq.T[tmp_x, tmp_y1:tmp_y2] = ndimage.shift(img_dq.T[tmp_x, tmp_y1:tmp_y2], -(y_on_G_img[j]%1), order = 1, mode ='wrap')
tmp_spec_2d[i] = img.T[tmp_x][tmp_y1+1:tmp_y2-1]
tmp_wht_2d[i] = img_wht.T[tmp_x][tmp_y1+1:tmp_y2-1]
tmp_dq_2d[i] = img_dq.T[tmp_x][tmp_y1+1:tmp_y2-1]
elif pupil == 'C':
tmp_y, tmp_x1, tmp_x2 = int(y_on_G_img[j]), int(x_on_G_img[j] - aper_int_pix - 1), int(x_on_G_img[j] + aper_int_pix + 2)
img[tmp_y, tmp_x1:tmp_x2] = ndimage.shift(img[tmp_y, tmp_x1:tmp_x2], -(x_on_G_img[j]%1), order = 1, mode ='wrap')
img_wht[tmp_y, tmp_x1:tmp_x2] = ndimage.shift(img_wht[tmp_y, tmp_x1:tmp_x2], -(x_on_G_img[j]%1), order = 1, mode ='wrap')
img_dq[tmp_y, tmp_x1:tmp_x2] = ndimage.shift(img_dq[tmp_y, tmp_x1:tmp_x2], -(x_on_G_img[j]%1), order = 1, mode ='wrap')
tmp_spec_2d[i] = img[tmp_y, tmp_x1+1:tmp_x2-1]
tmp_wht_2d[i] = img_wht[tmp_y, tmp_x1+1:tmp_x2-1]
tmp_dq_2d[i] = img_dq[tmp_y, tmp_x1+1:tmp_x2-1]
tmp_spec_2d = tmp_spec_2d.T
tmp_wht_2d = tmp_wht_2d.T
tmp_dq_2d = tmp_dq_2d.T
'''Construct fits file for 2d grism spectra'''
### Primary HDU
hdu = fits.PrimaryHDU()
hdu.header['x0'] = (np.float32(x0), 'Reference position X in direct image')
hdu.header['y0'] = (np.float32(y0), 'Reference position Y in direct image')
# tmp_coord = wcs.utils.pixel_to_skycoord(x0, y0, tmp_wcs)
# hdu.header['RA0'] = (tmp_coord.ra.value, 'Reference position RA in direct image')
# hdu.header['DEC0'] = (tmp_coord.dec.value, 'Reference position Dec in direct image')
hdu.header['author'] = ('Fengwu Sun', 'Author of this file')
hdu.header['time'] = (time.strftime("%Y/%m/%d %H:%M:%S", time.localtime()), 'Time of Creation')
### Scientific HDU
hdu_sci = fits.ImageHDU(np.float32(tmp_spec_2d), name = 'SPEC2D')
hdu_sci.header['wave_1'] = (wave[args_eff[0][1]], 'Wavelength (um) of first pixel')
hdu_sci.header['d_wave'] = (np.mean(np.diff(wave[args_eff[0]])), 'Wavelength Difference (um) between each pixel')
hdu_sci.header['comments'] = ('wave = wave_1 + np.arange(0, NAXIS1, 1) * d_wave')
# hdu_sci.header['pixscale'] = (pix_scale, 'Pixel scale in undispersed direction (arcsec)')
hdu_sci.header['aperture'] = (aper, 'Aperture radius in undispersed direction (arcsec)')
hdu_sci.header['pupil'] = (pupil, 'Pupil of grism (R or C)')
hdu_sci.header['module'] = (header['module'], 'Module of Detector (A or B)')
hdu_sci.header['diff_y'] = (y_on_G_img[args_eff[0]][0], 'Y_(full)-Y_(trim)')
hdu_sci.header['diff_x'] = (x_on_G_img[args_eff[0]][0], 'X_(full)-X_(trim)')
### Weight / Data-quality HDU
hdu_wht = fits.ImageHDU(np.float32(tmp_wht_2d), name = 'WHT2D')
hdu_dq = fits.ImageHDU(np.int32(tmp_dq_2d), name = 'DQ2D')
### Tracing & Dispersion information
tmp_tb_wave = Table(data = [wavs[args_eff], x_on_G_img[args_eff], y_on_G_img[args_eff], dxs[args_eff], dys[args_eff]],
names = ['wavelength', 'xs', 'ys', 'dxs', 'dys'])
tmp_tb_wave['wavelength'].info.format = '.6f'
for x in tmp_tb_wave.colnames[1:]: tmp_tb_wave[x].info.format = '.3f'
hdul = fits.HDUList([hdu, hdu_sci, hdu_wht, hdu_dq])
hdul.append(fits.BinTableHDU(tmp_tb_wave, name = 'WAVE'))
return hdul
def store_all_2d_spec(fits_list, pupils, modules, paths, output = 'test.fits',
coord = None, filter = None, info_table = None, overwrite = True):
'''
Write all the individual 2D spectra of one single sources into a fits file
-----------------------------------------------
Parameters
----------
fits_list : list of `astropy.io.fits`
list of fits files generated by extract_2d_spec()
pupils: list of strings
pupil list of grism ('R' or 'C')
modules: list of strings
module list of NIRCam detector ('A' or 'B')
paths: list of strings
paths of grism exposures
# poms: list of integers
# list of POM values
output: string
Output file name
coord: `astropy.coordinates.SkyCoord` object
coordinates of the source. default: None (will not record this in the header)
filter: string
Filter of grism observation (e.g., 'F444W', 'F332W', 'F356W')
default: None (will not record this in the header)
info_table: astropy.table.Row
A row of information related to this source
default: None (will not record this in the header)
overwrite: bool
If true, save (and overwrite) the data. default: True
Returns
-------
ind_hdul: `~astropy.io.fits.HDUList`
fits HDU list of individual sci, wht and dq data.
'''
## HDU list of individual exposures:
ind_hdul = fits.HDUList([fits.PrimaryHDU()])
for l, x in enumerate(fits_list):
ind_hdul.append(x[1])
ind_hdul[-1].header['EXTNAME'] = 'SPEC2D-%d' % l
for card_name in ['x0', 'y0']: # 'ra0', 'dec0']:
ind_hdul[-1].header[card_name] = x[0].header[card_name]
ind_hdul[-1].header['pupil'] = pupils[l]
ind_hdul[-1].header['module'] = modules[l]
# ind_hdul[-1].header['pom'] = poms[l]
ind_hdul[-1].header['datapath'] = paths[l].split('/')[-1]
ind_hdul.append(x[2])
ind_hdul[-1].header['EXTNAME'] = 'WHT2D-%d' % l
ind_hdul.append(x[3])
ind_hdul[-1].header['EXTNAME'] = 'DQ2D-%d' % l
ind_hdul.append(x[4])
ind_hdul[-1].header['EXTNAME'] = 'WAVE-%d' % l
#### stats of single exposures
ind_hdul.append(fits.BinTableHDU(Table(names = ['id', 'pupil', 'module', 'datapath'], # 'POM',
data = [range(len(fits_list)), pupils, modules, # poms,
[tmp_path.split('/')[-1] for tmp_path in paths]]),
name = 'STATS'))
#### add other important specs
ind_hdul[0].header['DIRNAME'] = (os.path.dirname(paths[0]), 'Directory name of original grism data')
if coord != None:
ind_hdul[0].header['RA0'] = (tmp_coord.ra.value, 'Reference position RA in direct image')
ind_hdul[0].header['DEC0'] = (tmp_coord.dec.value, 'Reference position Dec in direct image')
ind_hdul[0].header['N_coadd'] = (len(ind_hdul[-1].data['pupil']), 'Number of coadded frames')
ind_hdul[0].header['N_R'] = (np.sum(ind_hdul[-1].data['pupil'] == 'R'), 'Number of R grism frames')
ind_hdul[0].header['N_C'] = (np.sum(ind_hdul[-1].data['pupil'] == 'C'), 'Number of C grism frames')
ind_hdul[0].header['author'] = ('Fengwu Sun', 'Author of this file')
ind_hdul[0].header['time'] = (time.strftime("%Y/%m/%d %H:%M:%S", time.localtime()), 'Time of Creation')
if filter != None:
ind_hdul[0].header['filter'] = (tmp_filter, 'Filter name')
if info_table != None:
ind_hdul[0].header['COMMENTS'] = 'Belows are information taken from input catalog:'
for x in info_table.colnames:
ind_hdul[0].header['HIERARCH ' + x] = info_table[x]
if overwrite: ind_hdul.writeto(output, overwrite = overwrite) ## save the fits file?
return ind_hdul
def resample_spec2d_wmin_wmax(x, i):
'''
Resample a series of 2D spectra at wavelegth range = tmp_w_min - tmp_w_max
-----------------------------------------------
Parameters
----------
x: int
index of the individual fits files in the list.
i: int
index of the resampled wave array
tmp_wave_2d:
Global variable (`list(np.array())`).
A compilation of 1D wavelength array associated with all individual grism frames.
tmp_ind_fits_list:
Global variable (`astropy.io.fits.HDUList()`).
A HDU List that have extracted 2D spectra (SCI/WHT/DQ) from each individual frames
wave_sample:
Global variable (`np.array()`)
Resampled wavelengths.
Returns
-------
tmp_spec_w_1: `~np.array()`
a column of 2D SCI spectra from frame [x] that fall in the resampled wavelength range [i]
tmp_wht_w_1: `~np.array()`
a column of 2D WHT spectra at the same position as that of `tmp_spec_w_1`
tmp_cov_w_1: `~np.array()`
a column of 2D DQ spectra at the same position as that of `tmp_spec_w_1`
'''
global tmp_wave_2d, tmp_ind_fits_list, wave_sample # arrays to use
# global arr_spec_2d, arr_wht_2d, arr_cov_2d # arrays to save; impossible with pool
tmp_w_min, tmp_w_max = wave_sample[i], wave_sample[i+1]
## resampled SCI, WHT, DQ image
arg_in_wmin_wmax = tuple([(tmp_wave_2d[x] > tmp_w_min) & (tmp_wave_2d[x] <= tmp_w_max)])
tmp_spec_w = tmp_ind_fits_list[x][1].data.T[arg_in_wmin_wmax]
tmp_wht_w = np.nan_to_num(tmp_ind_fits_list[x][2].data.T[arg_in_wmin_wmax], posinf = 0, neginf = 0)
tmp_dq_w = tmp_ind_fits_list[x][3].data.T[arg_in_wmin_wmax]
tmp_wht_w[tmp_dq_w!=0] = 0
## stack/resample N spectra at wmin-wmax (M rows)
tmp_spec_w_1 = np.nansum(tmp_spec_w * tmp_wht_w, axis = 0) / np.nansum(tmp_wht_w, axis = 0) #/ sum_tmp_wht_w #
tmp_wht_w_1 = np.nansum(tmp_wht_w, axis = 0)
tmp_cov_w_1 = np.int8(tmp_wht_w_1 != 0)
## save resampled spectra to array
# arr_spec_2d[i, x], arr_wht_2d[i, x], arr_cov_2d[i, x] = tmp_spec_w_1, tmp_wht_w_1, tmp_cov_w_1
return(tmp_spec_w_1, tmp_wht_w_1, tmp_cov_w_1)
In [ ]:
Extract 2D grism spectra with the loop below:¶
In [68]:
'''Choose extraction mode:
all : R + C combined
R : R only
C : C only
'''
extract_mode = 'all' # all, R, C
start = time.time()
for i, idx in enumerate(Index_interest[:]):
'''
Read source information from the input catalog
'''
t_a = time.time()
## item_sex = tb_sex[tb_sex['NUMBER'] == idx][0]
## tmp_mag = item_sex['MAG_AUTO']
item_sex = tb_sex[tb_sex['id'] == idx][0]
tmp_mag = item_sex['abmag']
tmp_coord = SkyCoord(item_sex['RA'], item_sex['DEC'], unit = (u.deg, u.deg))
### tmp_coord = SkyCoord(item_sex['X_WORLD'], item_sex['Y_WORLD'], unit = (u.deg, u.deg))
print('-[%3d] ID%-5d mag=%.2f' % (i, idx, tmp_mag), end = " ")
'''
### if no record of grism exposure is found in the dict `grism_frame_ID_per_source`,
### the script will continue to save time
'''
if len(grism_frame_ID_per_source[str(idx)]) == 0: print(' >> no spec found; '); continue
else:
## find 2D grism image, get the list
this_list_rate_wcs = np.array(grism_frame_path_per_source['path_'+str(idx)])
this_ID_rate_wcs = np.array(grism_frame_ID_per_source[str(idx)])
# print('[%d spectra]' % len(this_list_rate_wcs))
### but, if the number of available frame is too small, we will skip the extraction
if len(grism_frame_ID_per_source[str(idx)]) <= 2: print(' >> no spec found; '); continue
'''
Examine each reduced image to extract 2D spectra
'''
tmp_ind_fits_list = []
tmp_ind_fits_name_list = []
tmp_ind_module_list = []
tmp_ind_pupil_list = []
tmp_ind_pom_list = []
for j, tmp_grism in enumerate(this_list_rate_wcs):
### `tmp_grism` is the path of grism data ('*/jw*.fits')
# if idx not in list_sl_table[j]['Index'] : continue ## legacy command
### `tmp_idx_grism` is the sequential number of grism data (like 0-95 with 96 frames)
tmp_idx_grism = this_ID_rate_wcs[j]
tmp_grism_img = fits.getdata(tmp_grism, 'sci') ### Grism SCI image data
tmp_grism_dq = fits.getdata(tmp_grism, 'dq') ### Grism DQ data
tmp_grism_hd = fits.getheader(tmp_grism) ### Grism primary Header
tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]
### depending on `extract_mode`, only extrac R/C data
if extract_mode.lower() == 'r':
if tmp_pupil == 'C' : continue
if extract_mode.lower() == 'c':
if tmp_pupil == 'R' : continue
tmp_grism_hd_sci = fits.getheader(tmp_grism, 'sci') ### Grism SCI Header
### Read or Generate Grism WHT data:
tmp_grism_wht_path = tmp_grism.replace('lv1.5.fits', 'wht.fits')
if os.path.isfile(tmp_grism_wht_path):
tmp_grism_wht = fits.getdata(tmp_grism_wht_path)
else:
### 2D grism WHT file not found, generate from Error map:
tmp_grism_err = fits.getdata(tmp_grism, 'err')
tmp_grism_err[tmp_grism_err == 0] = np.nan
tmp_grism_wht = tmp_grism_err**-2
tmp_grism_wht_hdul = fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(tmp_grism_wht, name = 'WHT')])
tmp_grism_wht_hdul[1].header = tmp_grism_hd_sci
tmp_grism_wht_hdul.writeto(tmp_grism_wht_path, overwrite = True)
tmp_grism_wht_hdul.close()
### Direct imaging positions
item_sl = list_sl_table[tmp_idx_grism][list_sl_table[tmp_idx_grism]['Index'] == idx][0]
x0 = item_sl['pixel_x']
y0 = item_sl['pixel_y']
print('%s(%.1f, %.1f) ' % (tmp_pupil, x0, y0), end = ' ')
# if (x0 > 2100) & (idx == 1357) & (tmp_pupil == 'R'): continue
### spectral tracing parameters
fit_opt_fit = list_fit_opt_fit[list_mod_pupil == tmp_module + tmp_pupil][0]
w_opt = list_w_opt[list_mod_pupil == tmp_module + tmp_pupil][0]
dxs, dys, wavs = grism_conf_preparation(x0 = x0, y0 = y0, pupil = tmp_pupil,
fit_opt_fit = fit_opt_fit, w_opt = w_opt)
tmp_aper = 15. # aperture size
### Extract 2D spectrum on each single frame
try:
tmp_ind_fits = extract_2d_spec(img = tmp_grism_img, WRANGE = WRANGE,
x0 = x0, y0 = y0, dxs = dxs, dys = dys,
wave = wavs, img_wht = tmp_grism_wht, img_dq = tmp_grism_dq,
header = tmp_grism_hd, pupil = tmp_pupil, aper = tmp_aper)
except ValueError:
continue
tmp_ind_fits_list.append(tmp_ind_fits)
tmp_ind_module_list.append(tmp_module)
tmp_ind_fits_name_list.append(tmp_grism)
tmp_ind_pupil_list.append(tmp_pupil)
### If no spectra can be extract, continue; else: save spectra of all individal extraction
if len(tmp_ind_fits_list) == 0:
print(' >> no spec found; ')
continue
else:
print(' >> N=%d, x0=%.1f, y0=%.1f' % (len(tmp_ind_fits_list), x0, y0), end = '')
if extract_mode.lower() not in ['r', 'c']:
'''Change the path of saved 2D spectra: '''
tmp_ind_specs_name = '/xdisk/egami/fengwusun/ERSs/IceAge/extract_2d/allspec_2d_%s_%s_ID%s.fits' % (tmp_target, tmp_filter, idx)
store_all_2d_spec(fits_list = tmp_ind_fits_list, pupils = tmp_ind_pupil_list,
modules = tmp_ind_module_list, paths = tmp_ind_fits_name_list,
output = tmp_ind_specs_name, coord = tmp_coord,
filter = tmp_filter, info_table = item_sex, overwrite = True)
print(' >> save all extracted spec2d ', end = '')
### Correct for sensitivity:
for k, tmp_fits in enumerate(tmp_ind_fits_list[:]):
tmp_modpup = tmp_fits[1].header['MODULE'] + tmp_fits[1].header['PUPIL']
tmp_f_sens = list_f_sens[np.where(list_mod_pupil == tmp_modpup)[0][0]]
tmp_wavelength = tmp_fits[4].data['wavelength']
tmp_ind_fits_list[k][1].header['bunit'] = ('mJy', 'Brightness Unit')
tmp_ind_fits_list[k][1].data = tmp_fits[1].data / tmp_f_sens(tmp_wavelength) * 1e3 # to unit of mJy
tmp_ind_fits_list[k][2].data = tmp_fits[2].data * tmp_f_sens(tmp_wavelength)**2 * 1e-6
'''
Combine 2D spectra from all extraction
'''
tmp_wave_2d = np.array([x[-1].data['wavelength'] for x in tmp_ind_fits_list], dtype = object)
### resampled to the following wavelength frame:
d_wave = 0.001
wave_sample = np.arange(WRANGE[0], WRANGE[1] + d_wave, d_wave)
wave_sample_c = wave_sample[:-1] + d_wave / 2.
arr_spec_2d = np.zeros((len(wave_sample)-1, len(tmp_wave_2d), len(tmp_ind_fits_list[0][1].data)))
arr_wht_2d = np.zeros_like(arr_spec_2d)
arr_cov_2d = np.zeros_like(arr_spec_2d)
np.seterr(invalid = 'ignore')
### New method: resample spectrum using a multiple-processing tool:
xx, ii = np.meshgrid(range(len(tmp_wave_2d)), range(len(wave_sample) - 1))
xx = xx.flatten()
ii = ii.flatten()
with Pool(np.min([8, len(xx)])) as pool:
output = pool.starmap(resample_spec2d_wmin_wmax, zip(xx, ii))
for o in range(len(output)):
i, x = ii[o], xx[o]
arr_spec_2d[i, x], arr_wht_2d[i, x], arr_cov_2d[i, x] = output[o]
## sigma clip spectra SED
tmp_sigma_clip = SigmaClip(sigma_lower = 2.0, sigma_upper = 3.0)
arg_sigclip = np.where(tmp_sigma_clip(np.nanmean(arr_spec_2d, axis = (0, 2))).mask == False)[0]
tmp_ind_fits_list = [tmp_ind_fits_list[x] for x in arg_sigclip]
tmp_ind_module_list = np.array(tmp_ind_module_list)[arg_sigclip]
tmp_ind_fits_name_list = np.array(tmp_ind_fits_name_list)[arg_sigclip]
arr_spec_2d = arr_spec_2d[:,arg_sigclip,:]
arr_wht_2d = arr_wht_2d[:,arg_sigclip,:]
arr_cov_2d = arr_cov_2d[:,arg_sigclip,:]
## weighted mean 2d spectra:
tmp_spec_2d = np.nansum(arr_spec_2d * arr_wht_2d, axis = 1) / np.nansum(arr_wht_2d, axis = 1) # sum_arr_wht_2d #
tmp_wht_2d = np.nansum(arr_wht_2d, axis = 1)
tmp_cov_2d = np.nansum(arr_cov_2d, axis = 1)
'''
Save Coadded 2D spectra
'''
tmp_tb_cov = Table(names = ['index', 'name', 'x0', 'y0', 'module', 'pupil', 'DIFF_X', 'DIFF_Y', 'wave_0', 'wave_1'],
data = [range(len(tmp_ind_fits_list)),
[tmp_grism.split('/')[-1][:-16] for tmp_grism in tmp_ind_fits_name_list],
[x[0].header['x0'] for x in tmp_ind_fits_list],
[x[0].header['y0'] for x in tmp_ind_fits_list],
tmp_ind_module_list,
[x[1].header['PUPIL'] for x in tmp_ind_fits_list],
[x[1].header['DIFF_X'] for x in tmp_ind_fits_list],
[x[1].header['DIFF_Y'] for x in tmp_ind_fits_list],
[np.round(wave_sample_c[np.sum(arr_wht_2d, axis = -1)[:, n_] != 0][0], 4)
for n_ in range(len(tmp_ind_fits_list))],
[np.round(wave_sample_c[np.sum(arr_wht_2d, axis = -1)[:, n_] != 0][-1], 4)
for n_ in range(len(tmp_ind_fits_list))],
])
hdu = fits.PrimaryHDU()
hdu.header['RA0'] = (tmp_coord.ra.value, 'Reference position RA in direct image')
hdu.header['DEC0'] = (tmp_coord.dec.value, 'Reference position Dec in direct image')
hdu.header['N_coadd'] = (len(tmp_tb_cov['pupil']), 'Number of coadded frames')
hdu.header['N_R'] = (np.sum(tmp_tb_cov['pupil'] == 'R'), 'Number of R grism frames')
hdu.header['N_C'] = (np.sum(tmp_tb_cov['pupil'] == 'C'), 'Number of C grism frames')
hdu.header['author'] = ('Fengwu Sun', 'Author of this file')
hdu.header['time'] = (time.strftime("%Y/%m/%d %H:%M:%S", time.localtime()), 'Time of Creation')
hdu.header['filter'] = (tmp_filter, 'Filter name')
# hdu.header['abmag'] = (tmp_mag, 'AB Magnitude in %s Band' % tmp_filter)
hdu.header['COMMENTS'] = 'Belows are information taken from input mock catalog:'
for x in item_sex.colnames:
hdu.header['HIERARCH ' + x] = item_sex[x]
hdu_sci = fits.ImageHDU(tmp_spec_2d.T, name = 'SPEC2D')
hdu_sci.header['wave_1'] = (wave_sample_c[0], 'Wavelength (um) of first pixel')
hdu_sci.header['d_wave'] = (wave_sample_c[1] - wave_sample_c[0], 'Wavelength Difference (um) between each pixel')
hdu_sci.header['comments'] = ('wave = wave_1 + np.arange(0, NAXIS1, 1) * d_wave')
hdu_sci.header['pixscale'] = (0.0629, 'Pixel scale in undispersed direction (arcsec)')
hdu_sci.header['aperture'] = (tmp_aper, 'Aperture radius in undispersed direction (arcsec)')
hdu_wht = fits.ImageHDU(tmp_wht_2d.T, name = 'WHT2D')
hdu_wht.header['comments'] = ('Weight image; ERR = WHT^(-0.5)')
hdu_cov = fits.ImageHDU(tmp_cov_2d.T, name = 'COV2D')
hdu_cov.header['comments'] = ('Coverage image')
hdu_tab = fits.BinTableHDU(tmp_tb_cov, name = 'STATS')
hdu_tab.header['COMMENT'] = 'name: name of simulated image'
hdu_tab.header['COMMENT'] = 'x0 / y0: reference position (i.e., in direct image)'
hdu_tab.header['COMMENT'] = ' (in reality this is registered to the wcs of grism image)'
hdu_tab.header['COMMENT'] = 'DIFF_X: X_(full)-X_(trim)'
hdu_tab.header['COMMENT'] = 'DIFF_Y: Y_(full)-Y_(trim)'
hdu_tab.header['COMMENT'] = 'wave_0: minimum wavelength (micron) in this coverage'
hdu_tab.header['COMMENT'] = 'wave_1: maximum wavelength (micron) in this coverage'
hdul = fits.HDUList([hdu, hdu_sci, hdu_wht, hdu_cov, hdu_tab])
'''Change the path of saved 2D spectra: '''
tmp_spec_2d_name = '/xdisk/egami/fengwusun/ERSs/IceAge/extract_2d/spec_2d_%s_%s_ID%s_comb.fits' % (tmp_target, tmp_filter, idx)
if extract_mode.lower() in ['r', 'c']:
tmp_spec_2d_name = tmp_spec_2d_name.replace('_comb.fits', '_%s.fits' % extract_mode.upper())
hdul.writeto(tmp_spec_2d_name, overwrite = True)
print(' >> save stacked spec2d; ')
end = time.time()
print('run time = %.1fs = %.1fmin' % (end - start, (end - start)/60))
In [ ]:
In [ ]:
7. Examine 2D & 1D spectra in plots¶
With the 2D grism spectra you just extracted and combined above, it is now the time to examine the quality of extracted spectra. You can plot the 2D, 1D spectrum and direct imaging data in the same frame.
In [78]:
'''List of 2D co-added grism spectra that we just extracted:'''
list_spec = np.array(glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/extract_2d/spec_2d_CHAMMS1-CO_F444W_ID*_comb.fits'))
print(len(list_spec), 'spectra files in total.')
'''Direct imaging data
for IceAge data, it is avaialble at:
http://gxn.as.arizona.edu/~sunfengwu/IceAge/mosaics/IceAge_CHAMMS1-C2-FIELD_lw_F430M_visitall_modall_i2d.fits.tar.gz
'''
path_LW = '/xdisk/egami/fengwusun/ERSs/IceAge/all_mosaics/IceAge_CHAMMS1-C2-FIELD_lw_F430M_visitall_modall_i2d.fits'
img_LW = fits.getdata(path_LW, 'sci')
hd_LW = fits.getheader(path_LW, 'sci')
wcs_LW = wcs.WCS(hd_LW)
pix_LW = wcs.utils.proj_plane_pixel_area(wcs_LW) ** 0.5 * 3600.
In [ ]:
In [92]:
# coords_bound = SkyCoord(['11:06:00.7328 -77:21:29.267', '11:06:51.4474 -77:22:46.359'], unit = (u.h, u.deg))
# opt_bound = np.polyfit(coords_bound.ra.value, coords_bound.dec.value, deg = 1)
for i, tmp_spec_path in enumerate(list_spec[1:2]):
tmp_spec_fits = fits.open(tmp_spec_path)
'''source information '''
tmp_id = tmp_spec_fits[0].header['ID'] # tmp_spec_fits[0].header['NUMBER']
tmp_filter = tmp_spec_fits[0].header['FILTER'] # tmp_spec_fits[0].header['select'] #
tmp_mag = tmp_spec_fits[0].header['ABMAG'] # tmp_spec_fits[0].header['MAG_AUTO']
# if tmp_mag < 23.5: continue
tmp_flux_mJy = 10**(-0.4 * tmp_spec_fits[0].header['ABMAG']) * 3631e3 # 10**(-0.4 * tmp_spec_fits[0].header['MAG_AUTO']) * 3631e3
tmp_N_R = tmp_spec_fits[0].header['N_R']
tmp_N_C = tmp_spec_fits[0].header['N_C']
tmp_RA, tmp_DEC = tmp_spec_fits[0].header['RA0'], tmp_spec_fits[0].header['DEC0']
tmp_x0, tmp_y0 = tmp_spec_fits[0].header['xcentroid'] , tmp_spec_fits[0].header['ycentroid']
# tmp_x0, tmp_y0 = tmp_spec_fits[0].header['X_IMAGE'] - 1, tmp_spec_fits[0].header['Y_IMAGE'] - 1
# if np.polyval(opt_bound, tmp_RA) > tmp_DEC: tmp_module = 'A'
# else: tmp_module = 'B'
'''2d spec'''
tmp_spec_2d = tmp_spec_fits[1].data
tmp_wht_2d = tmp_spec_fits[2].data
wave_sample_c = tmp_spec_fits[1].header['WAVE_1'] + np.arange(tmp_spec_fits[1].header['NAXIS1']) * tmp_spec_fits[1].header['D_WAVE']
## if not enough wavelength coverage, then skip
if np.sum(np.isnan(np.nansum(tmp_spec_2d, axis = 0)) == False) < 200: continue
## get y center:
tmp_spec_ydir = np.nansum(tmp_spec_2d * tmp_wht_2d, axis = 1) / np.nansum(tmp_wht_2d, axis = 1)
tmp_yc = np.argmax(tmp_spec_ydir[12:19]) + 12
# tmp_yc = 15
tmp_aper = 3
## model background with 1d fit at each x-pix position
h_space_bkg = np.append(np.arange(0, tmp_yc-6), np.arange(tmp_yc+6, 31))
tmp_specbkg_2d = tmp_spec_2d[h_space_bkg]
_, tmp_med_specbkg_2d, tmp_rms_specbkg_2d = sigma_clipped_stats(tmp_specbkg_2d, axis = 0)
tmp_specbkg_2d_3sigma = np.ones_like(tmp_specbkg_2d) * (tmp_med_specbkg_2d + tmp_rms_specbkg_2d * 2.5)
tmp_wht_specbkg_2d = np.clip(np.float32(tmp_specbkg_2d < tmp_specbkg_2d_3sigma), 0.01, 1)
tmp_specbkg_2d_model = np.vstack([np.polyval(np.polyfit(h_space_bkg, np.nan_to_num(tmp_specbkg_2d[:, x]), deg = 1,
w = tmp_wht_specbkg_2d[:, x]),
np.arange(31)) for x in range(tmp_specbkg_2d.shape[1])])
tmp_specbkg_2d_model = tmp_specbkg_2d_model.T
tmp_spec_2d = tmp_spec_2d - np.nan_to_num(tmp_specbkg_2d_model)
tmp_spec_1d = np.nansum(tmp_spec_2d[tmp_yc-tmp_aper:tmp_yc+tmp_aper+1], axis = 0)
tmp_unc_1d = np.nansum(tmp_wht_2d[tmp_yc-tmp_aper:tmp_yc+tmp_aper]**-1, axis = 0)**0.5
tmp_unc_1d = ((tmp_rms_specbkg_2d / 4.)**2 + tmp_unc_1d**2)**0.5
'''
Plot Spectra:
'''
plt.close()
fig, ax = plt.subplots(2, 1, figsize = (12, 5))
### ax[0]: 2D spectra
ax[0].imshow(tmp_spec_fits[1].data[:], aspect = 5, vmax = np.clip(tmp_flux_mJy / 2., 0.01, 1e5), vmin = -0.01,
cmap = plt.cm.gist_gray_r)
ax[0].set(ylim = (0.5, 30.5), xticks = [], aspect = 6 * (np.diff(WRANGE) - 0.1) / 1.2,
xlim = (0.05 / tmp_spec_fits[1].header['D_WAVE'],
(WRANGE[1] - WRANGE[0] - 0.05) / tmp_spec_fits[1].header['D_WAVE'])
)
ax[0].set_yticks([15.5])
ax[0].set_yticklabels([""])
### ax[1]: 1D spectra
ax[1].plot(wave_sample_c, ndimage.gaussian_filter1d(tmp_spec_1d, 0.5), # tmp_spec_1d, #
color = 'k', zorder = 100, lw = 1.5)
## uncomment below to add R/C-only spectra, if you have extracted them:
# tmp_spec_R_path = tmp_spec_path.replace('_comb.fits', '_R.fits')
# tmp_spec_C_path = tmp_spec_path.replace('_comb.fits', '_C.fits')
# if os.path.isfile(tmp_spec_R_path):
# tmp_r_fits = fits.open(tmp_spec_R_path)
# wave_sample_gR = tmp_r_fits[1].header['WAVE_1'] + np.arange(tmp_r_fits[1].header['NAXIS1']) * tmp_r_fits[1].header['D_WAVE']
# tmp_spec_R_1d = np.sum(tmp_r_fits[1].data[tmp_yc-3:tmp_yc+4], axis = 0)
# tmp_specbkg_R_2d = np.vstack((tmp_r_fits[1].data[0:10], tmp_r_fits[1].data[20:]))
# tmp_bkg_R_1d = sigma_clipped_stats(tmp_specbkg_R_2d, axis = 0, sigma_upper = 2)[0]
# tmp_spec_R_1d = tmp_spec_R_1d - tmp_bkg_R_1d * 10
# ax[1].plot(wave_sample_gR, tmp_spec_R_1d , color = 'r', zorder = 50, lw = 0.5)
# if os.path.isfile(tmp_spec_C_path):
# tmp_c_fits = fits.open(tmp_spec_C_path)
# wave_sample_gC = tmp_c_fits[1].header['WAVE_1'] + np.arange(0, tmp_c_fits[1].header['NAXIS1'], 1) * tmp_c_fits[1].header['D_WAVE']
# tmp_spec_C_1d = np.sum(tmp_c_fits[1].data[tmp_yc-3:tmp_yc+4], axis = 0)
# tmp_specbkg_C_2d = np.vstack((tmp_c_fits[1].data[0:10], tmp_c_fits[1].data[20:]))
# tmp_bkg_C_1d = sigma_clipped_stats(tmp_specbkg_C_2d, axis = 0, sigma_upper = 2)[0]
# tmp_spec_C_1d = tmp_spec_C_1d - tmp_bkg_C_1d * 10
# ax[1].plot(wave_sample_gC, tmp_spec_C_1d , color = 'c', zorder = 49, lw = 0.8)
tmp_max_counts = np.nanpercentile(tmp_spec_1d[(np.isnan(tmp_spec_1d) == False) & (tmp_spec_1d != 0)], 95) * 1.5
if np.isnan(tmp_max_counts): continue
ax[1].axhline(0, color = 'grey', ls = '--' )
ax[1].set(xlim = (WRANGE[0] + 0.05, WRANGE[1] - 0.05),
xticks = np.arange(WRANGE[0] + 0.1, WRANGE[1] - 0.05 + 0.01, 0.1),
xlabel = 'Observed Wavelength (µm)', ylabel = 'Flux Density [mJy]')
ax[1].set_ylim(-0.03, np.clip(tmp_max_counts, 0.03, 1e8))
corner_text(ax[0], loc = 2, s = 'ID%s' % (tmp_id), weight = 'semibold', color = 'r', fontsize = 18, edge = 5e-3)
corner_text(ax[0], loc = 1, s = '%s=%.2f' % ('F430M', tmp_mag), color = 'r', fontsize = 14, edge = 5e-3)
corner_text(ax[0], loc = 4, s = '(%.5f, %.5f)' % (tmp_RA, tmp_DEC), color = 'r', fontsize = 14, edge = 5e-3)
corner_text(ax[1], loc = 1, s = '$f_\mathrm{%s}$=%.3f mJy' % ('F430M', tmp_flux_mJy), fontsize = 14, color = 'r', edge = 5e-3)
corner_text(ax[1], loc = 2, s = 'N(R)=%2d' % tmp_N_R, fontsize = 12, color = 'r', edge = 5e-3)
corner_text(ax[1], loc = 3, s = 'N(C)=%2d' % tmp_N_C, fontsize = 12, color = 'c', edge = 5e-3)
fig.subplots_adjust(hspace = -0.20, wspace = 0, bottom = 0.13, top = 1.10, right = 0.98, left = 0.15)
### ax_im: Direct Image
ax_im = fig.add_axes([0.02, 0.68, 0.30 / 12 * 5, 0.30])
tmp_img = img_LW[int(tmp_y0)-16:int(tmp_y0)+17,int(tmp_x0)-16:int(tmp_x0)+17]
ax_im.imshow(tmp_img,
vmin = np.nanpercentile(tmp_img, 2.5), vmax = np.nanpercentile(tmp_img, 97.5),
cmap = plt.cm.gist_heat)
ax_im.set(xlim = (tmp_x0%1 + 0.5, tmp_x0%1 + 30.5), ylim = (tmp_y0%1 + 0.5, tmp_y0%1 + 30.5))
ax_im.set_xticks([])
ax_im.set_yticks([])
corner_text(ax_im, loc = 4, s = '(%.1f, %.1f)' % (tmp_x0, tmp_y0), color = 'w', fontsize = 10)
print(i, tmp_id, '%s=%.2f' % (tmp_filter, tmp_mag) )
tmp_spec_fits.close()
'''
Save spectral plot (remember to change the path)
'''
plt.savefig('/home/u24/fengwusun/IceAge/extracted_spectra/spec_1d_%s_%s_ID%s_mag%.2f.pdf' % (tmp_target, tmp_filter,
tmp_id, tmp_mag),
dpi = 150)
'''
Save 1D spectra (remember to change the path)
'''
tb_1dspec = Table(names = ['wavelength_um', 'flux_mJy', 'fluxerr_mJy'],
data = [wave_sample_c, tmp_spec_1d, tmp_unc_1d])
tb_1dspec['wavelength_um'].info.format = '.4f'
tb_1dspec['flux_mJy'].info.format = '.5f'
tb_1dspec['fluxerr_mJy'].info.format = '.5f'
ascii.write(tb_1dspec, tmp_spec_path.replace('/spec_2d_', '/spec_1d_').replace('.fits', '.dat'),
format = 'commented_header', overwrite = True)
Optional: Check individual source & multiple frames¶
You can overlie spectral tracing / dispersion models on different 2D grism frames & individual sources to check the accuracy the spectral tracing models, and apply shifts for the POM catalogs if necessary.
In [102]:
'''
Select frames for F444W grism images, plot tracing/dispersion models for top-20 bright sources:
'''
# glob.glob('/xdisk/egami/fengwusun/ERSs/IceAge/grism/jw0130902[5-6]00*_nrc[a-b]long_rate_lv1.5.fits')
for k, tmp_rate_path in enumerate(all_rate_list):
if k != 1: continue
### POM-applied catalog:
tmp_tb_sl = ascii.read('/home/u24/fengwusun/IceAge/POM_catalog/%s_%s_dirimg_sources.list' % (os.path.basename(tmp_rate_path)[:-20], tmp_filter))
### LW grism image:
tmp_rate_fits = fits.open(tmp_rate_path)
tmp_rate_wcs_base = tmp_rate_path.split('/')[-1].split('long_rate')[0]
tmp_grism_hd = tmp_rate_fits[0].header
tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]
print('[%-3d]' % k, tmp_rate_wcs_base)
print('/filter:', tmp_filter, '/module:', tmp_module, '/grism:', tmp_pupil)
### Source catalog - selected the brightest sources
tmp_tb_sl.sort('MAG_AUTO')
tmp_tb_sl = tmp_tb_sl[:20]
### Wave space for plotting
if tmp_filter == 'F322W2': wave_space = np.arange(2.4, 4.11, 0.2)
elif tmp_filter == 'F444W': wave_space = np.arange(3.85, 5.11, 0.2)
### Spectral tracing & Dispersion parameters:
if tmp_filter in ['F277W', 'F356W']: disp_filter = 'F322W2'
else: disp_filter = tmp_filter
tb_order23_fit = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, tmp_module, tmp_pupil))
fit_opt, fit_err = tb_order23_fit['col0'].data, tb_order23_fit['col1'].data
tb_fit_displ = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % (tmp_module, tmp_pupil))
w_opt, w_err = tb_fit_displ['col0'], tb_fit_displ['col1']
'''Plot 2D grism spectra'''
plt.close()
fig, ax = plt.subplots(1, 1, figsize = (12, 12))
plt.imshow(tmp_rate_fits[1].data, vmin = -0.03, vmax = np.percentile(tmp_rate_fits[1].data, 98.),
origin = 'lower', aspect = 1, cmap = plt.cm.viridis)
### show stars
for i, item_sl in enumerate(tmp_tb_sl):
if item_sl['Index'] == 1: kwargs_target = dict(marker = '*', color = 'pink', mec = 'w', mew = 2, ms = 20)
else:
kwargs_target = dict(marker = '*', color = 'r', mec = 'pink', mew = 2, ms = 20)
plt.plot(item_sl['pixel_x'], item_sl['pixel_y'], ls = 'none', **kwargs_target)
dx_space = func_fit_wave(np.vstack((item_sl['pixel_x'] * np.ones_like(wave_space),
item_sl['pixel_y'] * np.ones_like(wave_space), wave_space)), *w_opt)
dy_space = fit_disp_order23(np.vstack((item_sl['pixel_x'] * np.ones_like(dx_space),
item_sl['pixel_y'] * np.ones_like(dx_space), dx_space)), *fit_opt)
if tmp_pupil == 'R':
plt.plot(item_sl['pixel_x'] + dx_space, item_sl['pixel_y'] + dy_space,
color = 'r', lw = 0.8, marker = '|', ms = 10, ls = ':')
for j, tmp_dx in enumerate(dx_space):
if ((item_sl['pixel_x'] + dx_space[j] > 2042) | (item_sl['pixel_x'] + dx_space[j] < 5)|
(item_sl['pixel_y'] + dy_space[j] > 2042) | (item_sl['pixel_y'] + dy_space[j] < 5)): continue
plt.text(item_sl['pixel_x'] + dx_space[j], item_sl['pixel_y'] + dy_space[j] - 3 ,
s = '%.1f' % wave_space[j], color = 'r', ha = 'center', va = 'top', fontsize = 10)
elif tmp_pupil == 'C':
plt.plot(item_sl['pixel_x'] + dy_space, item_sl['pixel_y'] + dx_space,
color = 'r', lw = 0.8, marker = '_', ms = 10, ls = ':')
for j, tmp_dx in enumerate(dx_space):
if ((item_sl['pixel_y'] + dx_space[j] > 2042) | (item_sl['pixel_y'] + dx_space[j] < 5)|
(item_sl['pixel_x'] + dy_space[j] > 2042) | (item_sl['pixel_x'] + dy_space[j] < 5)): continue
plt.text(item_sl['pixel_x'] + dy_space[j], item_sl['pixel_y'] + dx_space[j] - 3 ,
s = '%.2f' % wave_space[j], color = 'r', ha = 'right', va = 'center', fontsize = 10)
ax.set(xlim = (0, 2048), ylim = (0, 2048))
corner_text(ax, loc = 3, s = '%s mod%s (%s)' % (tmp_filter, tmp_module, tmp_pupil),
color = plt.cm.Reds(0.05), weight = 'bold', fontsize = 28, edge = 1e-2)
plt.tight_layout()
### Uncomment below if you want to save the plot; remember to change the path.
# fig.savefig('/home/u24/fengwusun/IceAge/tracing/%s_%s_%s_%s.pdf' % (tmp_rate_wcs_base, tmp_filter, tmp_module, tmp_pupil),
# dpi = 100)
In [ ]:
In [104]:
'''
Select one frame of F444W grism images, plot tracing/dispersion models for one "random" bright source:
'''
k = 1
tmp_rate_path = all_rate_list[k]
### POM-applied catalog
tmp_tb_sl = ascii.read('/home/u24/fengwusun/IceAge/POM_catalog/%s_%s_dirimg_sources.list' % (os.path.basename(tmp_rate_path)[:-20], tmp_filter))
### 2D grism image
tmp_rate_fits = fits.open(tmp_rate_path)
tmp_rate_wcs_base = tmp_rate_path.split('/')[-1].split('long_rate')[0]
tmp_grism_hd = tmp_rate_fits[0].header
tmp_filter, tmp_module, tmp_pupil = tmp_grism_hd['filter'], tmp_grism_hd['module'], tmp_grism_hd['pupil'][-1]
print('[%-3d]' % k, tmp_rate_wcs_base)
print('/filter:', tmp_filter, '/module:', tmp_module, '/grism:', tmp_pupil)
### Source catalog - sorted by brightness, selected ones with relatively complete spectral coverage
tmp_tb_sl.sort('MAG_AUTO')
tmp_tb_sl = tmp_tb_sl[tmp_tb_sl['MAG_AUTO'] < 20.5]
if tmp_filter == 'F322W2':
if tmp_pupil == 'C':
tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_y'] - 1600) < 500)][:]
elif (tmp_module, tmp_pupil) == ('A', 'R'):
tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_x'] - 1700) < 500)][:]
elif (tmp_module, tmp_pupil) == ('B', 'R'):
tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_x'] - 400) < 500)][:]
elif tmp_filter == 'F444W':
if tmp_pupil == 'C':
tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_y'] - 700) < 600)][:]
elif (tmp_module, tmp_pupil) == ('A', 'R'):
tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_x'] - 700) < 600)][:]
elif (tmp_module, tmp_pupil) == ('B', 'R'):
tmp_tb_sl = tmp_tb_sl[(np.abs(tmp_tb_sl['pixel_x'] - 1300) < 600)][:]
tmp_tb_sl = tmp_tb_sl[:50]
### random source, or give an ID
i_random = np.random.randint(len(tmp_tb_sl))
# i_random = 62
tmp_tb_sl = tmp_tb_sl[i_random:i_random+1]
### Wavelength space
if tmp_filter == 'F322W2': wave_space = np.arange(2.4, 4.11, 0.2)
elif tmp_filter == 'F444W': wave_space = np.arange(3.85, 5.11, 0.2)
### Spectral tracing and dispersion parameters:
if tmp_filter in ['F277W', 'F356W']: disp_filter = 'F322W2'
else: disp_filter = tmp_filter
tb_order23_fit = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISP_%s_mod%s_grism%s.dat' % (disp_filter, tmp_module, tmp_pupil))
fit_opt, fit_err = tb_order23_fit['col0'].data, tb_order23_fit['col1'].data
tb_fit_displ = ascii.read('/home/u24/fengwusun/Comissioning/FS_grism_config_v2_202208/DISPL_mod%s_grism%s.dat' % (tmp_module, tmp_pupil))
w_opt, w_err = tb_fit_displ['col0'], tb_fit_displ['col1']
'''Plot 2D spectra cutout from the 2d grism image'''
plt.close()
fig, ax = plt.subplots(1, 1, figsize = (12, 5))
if tmp_pupil == 'C':
plt.imshow(tmp_rate_fits[1].data.T, vmin = -0.03, vmax = np.percentile(tmp_rate_fits[1].data, 96.),
origin = 'lower', aspect = 1, cmap = plt.cm.viridis, interpolation = 'gaussian')
elif tmp_pupil == 'R':
if (tmp_filter, tmp_module) == ('F322W2', 'B'): tmp_vmax = np.percentile(tmp_rate_fits[1].data, 80.)
else: tmp_vmax = np.percentile(tmp_rate_fits[1].data, 96.)
plt.imshow(tmp_rate_fits[1].data, vmin = -0.03, vmax = tmp_vmax,
origin = 'lower', aspect = 1, cmap = plt.cm.viridis, interpolation = 'gaussian')
### show stars
for i, item_sl in enumerate(tmp_tb_sl):
if item_sl['Index'] == 1: kwargs_target = dict(marker = '*', color = 'pink', mec = 'w', mew = 2, ms = 20)
else:
kwargs_target = dict(marker = '*', color = 'r', mec = 'pink', mew = 2, ms = 20)
plt.plot(item_sl['pixel_x'], item_sl['pixel_y'], ls = 'none', **kwargs_target)
dx_space = func_fit_wave(np.vstack((item_sl['pixel_x'] * np.ones_like(wave_space),
item_sl['pixel_y'] * np.ones_like(wave_space), wave_space)), *w_opt)
dy_space = fit_disp_order23(np.vstack((item_sl['pixel_x'] * np.ones_like(dx_space),
item_sl['pixel_y'] * np.ones_like(dx_space), dx_space)), *fit_opt)
if tmp_pupil == 'R':
plt.plot(item_sl['pixel_x'] + dx_space, item_sl['pixel_y'] + dy_space,
color = 'r', lw = 0.8, marker = '|', ms = 10, ls = ':')
for j, tmp_dx in enumerate(dx_space):
if ((item_sl['pixel_x'] + dx_space[j] > 2042) | (item_sl['pixel_x'] + dx_space[j] < 5)|
(item_sl['pixel_y'] + dy_space[j] > 2042) | (item_sl['pixel_y'] + dy_space[j] < 5)): continue
plt.text(item_sl['pixel_x'] + dx_space[j], item_sl['pixel_y'] + dy_space[j] - 3 ,
s = '%.1f' % wave_space[j], color = 'r', ha = 'center', va = 'top', fontsize = 10)
elif tmp_pupil == 'C':
plt.plot(item_sl['pixel_y'] + dx_space, item_sl['pixel_x'] + dy_space, # + 2.0,
color = 'r', lw = 0.8, marker = '|', ms = 10, ls = ':')
for j, tmp_dx in enumerate(dx_space):
if ((item_sl['pixel_y'] + dx_space[j] > 2042) | (item_sl['pixel_y'] + dx_space[j] < 5)|
(item_sl['pixel_x'] + dy_space[j] > 2042) | (item_sl['pixel_x'] + dy_space[j] < 5)): continue
plt.text(item_sl['pixel_y'] + dx_space[j], item_sl['pixel_x'] + dy_space[j] - 3 ,
s = '%.2f' % wave_space[j], color = 'r', ha = 'right', va = 'center', fontsize = 10)
### xlim/ylim of the plot
if tmp_pupil == 'R':
ax.set(xlim = (item_sl['pixel_x'] + np.min(dx_space) - 20, item_sl['pixel_x'] + np.max(dx_space) + 20 ),
ylim = (item_sl['pixel_y'] + np.median(dy_space) - 7.5, item_sl['pixel_y'] + np.median(dy_space) + 7.5 ),
aspect = 'auto')
elif tmp_pupil == 'C':
ax.set(xlim = (item_sl['pixel_y'] + np.min(dx_space) - 20, item_sl['pixel_y'] + np.max(dx_space) + 20 ),
ylim = (item_sl['pixel_x'] + np.median(dy_space) - 7.5, item_sl['pixel_x'] + np.median(dy_space) + 7.5 ),
aspect = 'auto')
corner_text(ax, loc = 3, s = '%s mod%s (%s)' % (tmp_filter, tmp_module, tmp_pupil),
color = plt.cm.Reds(0.05), weight = 'bold', fontsize = 28, edge = 1e-2)
corner_text(ax, loc = 1, s = tmp_rate_wcs_base, color = 'w')
# plt.tight_layout()
Out[104]:
In [105]:
tmp_tb_sl
Out[105]:
In [ ]:
In [ ]: